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Executive  Summary 


This  report  describes  the  results  of  applying  a  new  numerical  simulation  technique 
to  the  problem  of  scattering  from  subsurface  bubble  clouds  generated  by  spilling  or 
breaking  waves.  To  the  extent  that  these  bubble  clouds  can  be  modeled  by  homogeneous 
voids  or  weakly  permeable  objects,  the  scattering  characteristics  are  amenable  to  exact 
analysis  by  our  method.  The  results  extend  previous  work  in  which  we  developed 
numerical  simulation  techniques  to  compute  detailed  scattering  characteristics  of  ocean- 
surface  realizations  that  preserve  the  essential  hydrodynamic  nonlinearities.  It  is  known 
that  at  high  sea  states  the  measured  surface  reverberation  noise  cannot  be  explained 
by  surface  scatter  alone,  and  subsurface  bubble  clouds  have  been  hypothesized  as  the 
source  of  the  enhanced  backscatter.  In  its  simplest  form,  the  model  invokes  a  meter- 
sized  impermeable  object,  which  can  have  a  low-frequency  acoustic  backscatter  cross 
section  that  liflnore  than  30  dB  higher  than  the  corresponding  surface  backscatter 
level  at  high  sea  states.  The  proximity  of  the  highly  irregular  surface  that  spawned 
the  bubble  cloud,  however,  can  significantly  alter  the  scattering  characteristics  of  the 
bubble  in  isolation. 

The  bubble  scattering  problem  is  solved  herein  as  a  special  case  of  a  more  general 
formulation  we  developed  for  calculating  the  multiple  scattering  from  a  collection  of 
objects  whose  free-space  scattering  characteristics  are  known.  Each  scatterer  is  con¬ 
fined  to  a  region  bounded  by  a  pair  of  planes  over  which  the  total  scattered  fields  are 
characterized  by  their  spatial  two-dimensional  Fourier  spectra.  At  each  boundary  plane 
these  wave  fields  see  only  the  two  adjacent  scatterers.  If  there  are  only  two  scatterers, 
the  interaction  equations  for  the  coupled  Fourier  modes  can  be  solved  as  a  system  of 
linear  equations.  The  technique  is  sufficiently  general  that  one  of  the  particles  can  be 
replaced  by  a  plane  or  rough  surface.  Thus,  the  mutual  interaction  between  an  object 
and  a  rough  surface  can  be  solved  in  terms  of  the  known  scattering  functions  for  the 
surface  and  the  particle.  For  the  ocean  surface,  we  use  our  previously  developed  nu¬ 
merical  method  to  calculate  the  surface  scattering  function.  At  the  present  time,  the 
bubble  cloud  is  modeled  by  a  simple  geometrical  shape  whose  scattering  characteristics 
are  known.  More  refined  bubble  scattering  models  can  be  incorporated  as  they  are 
developed. 

Our  simulations  are  presently  restricted  to  two  dimensions.  This  constraint  comes 
mainly  from  the  need  to  characterize  ocean  surface  scatter  at  high  sea  states.  No 
currently  available  computational  method  gives  the  two-dimensional  surface  scatter¬ 
ing  function  with  sufficient  accuracy  to  perform  the  mutual-interaction  computations. 
However,  both  cylinders  and  spheres  near  plane  impenetrable  surfaces  can  be  solved  by 
the  mutual-interaction  method.  Indeed,  for  small  objects  that  scatter  isotropically,  the 
mutual-interaction  equations  admit  algebraic  solutions.  We  find  that  as  a  small  scatter¬ 
ing  object  is  moved  away  from  a  plane  impenetrable  surface,  the  backscatter  coefficient. 
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which  is  initially  suppressed  to  the  surface  backscatter  value  if  the  roughness  is  small, 
increases  very  rapidly  to  its  free-space  value  and  then  exceeds  that  value  by  at  least  12 
dB  before  it  peaks.  Thus,  a  small  object  in  close  proximity  to  an  impenetrable  surface 
can  display  very  rapid  changes  in  its  overall  scatter  coefficient,  and  large  enhancements 
over  the  free  backscatter  coefficient  can  be  achieved.  The  backscatter  enhancement  is 
generally  maximized  for  isotropic  scatterers.  The  exact  solutions  for  meter-sized  objects 
show  peak  enhancements  2  to  7  dB  lower  than  the  isotropic  value  depending  on  the 
incidence  angle. 

While  the  plane  surface  and/or  isotropic  results  are  useful  for  illustrating  some  key 
features  of  the  scattering  interaction,  they  do  not  admit  a  direct  dependence  on  sea- 
state.  Acoustic  scatter  data  show  that  the  onset  of  anomalous  surface  backscatter  occurs 
at  a  wind  speed  of  approximately  20  knots.  Our  simulations  show  that  at  20  knots  the 
scattering  characteristics  are  qualitatively  similar  to  those  of  a  plane  reflecting  surface, 
but  at  the  depth  of  the  deepest  wave  trough  the  60°  backscatter  enhancement  is  already 
near  its  peak  level.  At  30  knots  there  are  significant  differences  between  the  flat  and 
rough-surface  models,  but  the  lobing  pattern  with  increasing  distance  persists.  As  the 
wind  speed  is  increased  further,  the  characteristic  interference  patterns  may  disappear 
altogether.  In  all  cases,  however,  at  high  sea  states  (at  or  above  15  knots),  the  presence 
of  a  meter-sized  void  acts  to  substantially  enhance  the  low-frequency  acoustic  surface 
reverberation  noise  at  the  level  of  the  the  deepest  wave  trough. 

The  mutual  interaction  method  presently  does  not  allow  the  bubble  cloud  and  the 
surface  to  overlap.  Thus,  the  bubble  cloud  must  be  at  or  below  the  depth  of  the 
deepest  trough  in  the  surface  realization.  Somewhat  fortuitously,  this  feature  provides 
a  qualitative  explanation  of  the  wind-speed  dependence  of  the  Chapman-Harris  curves. 
The  bubble  backscatter  coefficient  at  a  distance  corresponding  to  the  deepest  wave 
trough  has  the  correct  incidence-angle  dependence,  but  it  rises  too  steeply  for  a  fixed 
bubble  size  independent  of  wind  speed.  The  present  bubble  model  does  not  include 
dynamic  or  wind-speed  dependent  characteristics.  Thus,  we  have  not  computed  the 
Doppler  spectrum  because  the  static  model  will  automatically  concentrate  the  excess 
backscattered  intensity  near  zero  Doppler  shift.  It  is  more  appropriate  to  perform  the 
more  expensive  Doppler  computations  when  bubble  dynamics  can  be  included. 

The  main  purpose  of  this  report  is  to  provio  thorough  expose  of  the  mutual- 
interaction  method.  The  preliminary  results  do  no  •  establish  all  the  details  of  the 
backscatter  enhancement  caused  by  subsurface  bubble  clouds,  but  the  missing  features 
can  be  incorporated  as  the  bubble-cloud  dynamics  and  wind-speed  dependence  are  bet¬ 
ter  understood.  The  rationale  for  the  bubble  model  is  reviewed  in  Section  I.  Section  II 
presents  some  background  material,  including  cross  section  definitions.  In  Section  III.  1 
the  theory  is  developed  in  its  most  general  vector  form  and  then  specialized  to  the  scalar 
problem  of  interest  here.  Because  the  scattering  interactions  nre  exact,  it  is  necessary 
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to  accommodate  evanescent  waves.  Analytic  continuation  is  used  to  extend  the  usual 
definition  of  the  scattering  functions,  but  the  necessary  results,  which  are  not  readily 
available  in  papers  or  text  books,  are  developed  in  the  appendix.  The  formulation  that 
is  used  for  numerical  computation  is  summarized  in  Section  III. 2.  Illustrative  examples 
are  presented  in  Section  III. 3,  where  two-cylinder  results  are  compared  to  computations 
performed  with  full  method  of  moments.  Comparisons  are  also  made  with  analytic  for¬ 
mulas  valid  for  isotropic  scatterers. 

Our  preliminary  results  for  ocean  surfaces  are  presented  in  Section  IV.  Conclusions 
and  recommendations  are  presented  in  Section  V.  An  overview  of  the  report  can  be 
obtained  by  reading  Sections  I  and  V. 
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I  Introduction 

In  a  previous  report  [1]  we  described  numerical  simulations  of  acoustic  and  electromag¬ 
netic  (EM)  scatter  from  evolving  nonlinear  ocean-surface  realizations.  The  work  was 
undertaken  to  verify  and  extend  theories  that  attempt  to  explain  acoustic  scatter  from 
ocean  surfaces.  Of  particular  interest  is  an  anomalous  enhanced  backscatter  component 
observed  at  high  sea  states.  Because  ocean  surfaces  comprise  a  power-law  continuum 
of  spatial  wave  number  components,  neither  the  small  perturbation  theory  of  rough 
surfaces  nor  the  Kirchoff  approximation  can  be  used  to  predict  the  characteristics  of 
scattered  wave  fields.  The  two-scale  model  for  rough-surface  scattering  effectively  com¬ 
bines  these  two  approximations;  however,  it  is  not  entirely  satisfactory  because  the 
two-scale  partitioning  and  the  subsequent  analysis  of  the  scattering  regimes  are  ad  hoc. 
More  recent  theoretical  work  based  on  higher-order  perturbation  theory  [2]  and  a  T- 
matrix  approach  [3]  has  largely  removed  the  immediate  problems  with  the  two-scale 
model,  and  a  well-founded  theory  has  emerged  that  is  applicable  to  the  wavelengths 
and  sea  states  relevant  to  the  anomalous  acoustic  backscatter  phenomenon.  In  prin¬ 
ciple,  the  new  theory  can  accommodate  nonlinear  hydrodynamics  and  will  ultimately 
yield  average  Doppler  spectra,  but  this  has  not  yet  been  achieved. 

Our  previous  work  demonstrated  that  practical  schemes  exist  for  generating  sur¬ 
face  realizations  that  accommodate  the  dominant  nonlinear  surface  hydrodynamics. 
Numerical  simulations  that  use  these  realizations  as  inputs  are  effectively  computer 
experiments,  whereby  any  characteristic  of  the  scattered  signal,  such  as  its  Doppler 
spectrum,  is  readily  computed.  At  the  present  time,  however,  computation  and  stor¬ 
age  requirements  restrict  the  diffraction  computations  to  two  dimensions.  That  is,  the 
surface  height  can  vary  in  only  one  direction.  This  restricts  the  fidelity  of  the  surface 
hydrodynamic  and  scattering  phenomena  that  can  be  modeled,  but  the  two-dimensional 
simulations  can  be  used  effectively  to  test  hypotheses  and  to  guide  experiment  planning. 
Thus,  to  generate  timely  intermediate  results,  we  have  extended  our  two-dimensional 
simulations  to  accommodate  near-surface  bubble  clouds,  which  have  been  hypothesized 
to  explain  the  anomalous  acoustic  backscatter. 

From  first-order  perturbation  theory,  it  follows  that  the  acoustic  surface  backscatter 
level  is  proportional  to  the  spatial  wave  number  spectrum  of  the  surface- height  fluc¬ 
tuations  evaluated  at  the  Bragg  wave  number,  which  is  given  as  kb  =  2k  sin  where 
k  =  2n/X,  X  is  the  wavelength  of  the  acoustic  signal,  and  0,  is  the  incidence  angle.  For 
a  simple  spectral  model,  for  example,  the  Phillips  and  Moskowitz  spectrum  discussed  in 
Section  IV  of  [1],  the  Bragg  level  saturates  at  a  threshold  wind  speed.  The  more  detailed 
scatter  calculations  predict  a  nearly  linear  increase  in  the  logarithm  of  the  backscatter 
cross  section  with  wind  speed.  The  excess  backscatter  over  the  Bragg  level  is  attributed 
to  large-scale  tilts  upon  which  the  Bragg  carpet  lies.  Ocean-surface  scatter  measured 


1 


-20. 
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Figure  1:  Comparison  of  first-order  perturbation  theory  with  130-Hz  acoustic  backscat- 
ter  data. 

with  microwave  radars  is  generally  in  agreement  with  these  predictions,  but  at  high 
wind  speeds  the  acoustic  backscatter  cross  section  increases  more  nearly  as  the  fourth 
power  of  wind  speed.  When  this  anomalous  backscatter  occurs,  the  Doppler  spectra 
show  a  pronounced  enhancement  near  zero  Doppler  shift. 

Figure  1  shows  the  130-Hz  Bragg  acoustic  backscatter  level  at  70°  and  80°  incidence 
angles  as  predicted  by  first-order  perturbation  theory  with  the  Donelan- Pierson  spectral 
model  for  an  upwind  propagation  direction.  The  two  curves  marked  CH-70  and  CH- 
80  are  derived  from  analytic  approximations  to  the  Chapman-Harris  curves  [4],  which 
are  commonly  used  to  estimate  the  acoustic  surface  reverberation  noise.  As  noted 
above,  more  detailed  calculations  give  higher  surface  backscatter  levels,  but  for  winds 
above  approximately  20  knots,  surface  roughness  alone  cannot  account  for  the  observed 
scatter  levels.  Because  of  the  wind-speed  dependence,  it  has  been  hypothesized  that  the 
anomalous  backscatter  is  caused  by  subsurface  bubble  clouds  generated  by  breaking  or 
spilling  waves. 
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The  support  for  this  hypothesis  has  been  summarized  by  Ellinthorpe  [5].  From 
laboratory  measurements  by  Baldy  [6],  we  know  that  breaking  waves  generate  clouds  of 
small  bubbles  in  which  the  largest-size  constituents  extend  to  a  depth  of  approximately 
one  half  the  rms  wave  height  of  the  breaking  waves.  Below  this  depth  there  is  a  sharp 
cutoff,  and  only  very  small  bubbles  survive.  If  we  use  significant  wave  height  as  deduced 
from  the  Pierson-Neumann  theory  in  place  of  the  rms  wave  height,  the  cloud  of  larger 
bubbles  extends  to  the  depth 

dmnx  as  1.5  X  10 -2U2,  (1) 


where  U  is  the  wind  speed  in  meters  per  second  (cf.  Chapter  8.3  of  Kinsman  [7]).  Baldy 
also  observed  that  the  bubbles  tend  to  occur  in  clusters  or  clouds  at  one  per  dominant 
wave  period,  with  a  definite  position  relative  to  the  wave  phase.  If  we  use  the  peak  of 
the  Pierson-Moskowitz  spectrum  as  a  measure  of  the  dominant  wavelength,  the  bubble 
clouds  should  have  a  lateral  dimension  bounded  by 

*m»x  <  U2/(l0g),  (2) 


where  g  is  the  acceleration  due  to  gravity.  From  119-kHz  surface  acoustic  backscatter 
measurements  in  deep  ocean,  Crawford  and  Farmer  [8]  deduced  plume-like  clouds  of 
bubbles  that  support  the  general  features  of  Baldy’s  laboratory  measurements. 

The  structure  and  dynamics  of  subsurface  bubble  clouds  are  poorly  understood,  but 
we  can  deduce  likely  scattering  characteristics.  In  a  sparse  distribution  of  scatterers, 
the  coherent  wave  field  propagates  with  the  effective  wave  number 
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feeff  ~  k  +  —  </(ai,ai)>  , 


(3) 


where  /(a,,  aj  is  the  complex  scattering  function  for  the  particle,  and  the  angle  brackets 
denote  an  ensemble  average  over  the  particle  attributes-mainly  size  [9).  For  a  small 
bubble  the  scattering  is  nearly  isotropic,  and  it  is  shown  in  [9]  that 
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uj2  —  fi2  +  iwT  ’ 

where  a  is  the  bubble  radius,  uj  is  the  acoustic  angular  frequency, 


(4) 


(5) 


is  the  corresponding  resonant  frequency,  and  T  is  the  thermal  loss  factor.  The  parame¬ 
ters  p  and  c  denote  density  and  speed,  respectively,  with  the  subscripts  indicating  the 
medium  (air  or  water).  For  low  acoustic  frequencies  the  expected  bubble  radii  are  such 
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that  w«!l,  and  thermal  losses  are  negligible.  Indeed,  the  scattering  losses  as  inferred 
from  the  extinction  of  the  coherent  wave  field  are  themselves  very  small,  which  indi¬ 
cates  that  sea*' -ring  by  individual  bubbles  is  not  a  viable  mechanism  to  account  for 
anomalous  scatter. 

On  the  other  hand,  for  frequencies  well  below  it  can  be  shown  from  (3),  (4),  and 
(5)  that  the  effective  sound  speed  inside' the  bubble  cloud  is  given  by  the  relation 


css  0^(1  +  8000 VT), 


(6) 


where  cw  is  the  speed  of  sound  in  water  and  VF  is  the  void  fraction.  For  void  fractions 
in  the  expected  range  2  x  10~4  <  VF  <  2  x  10~3,  the  speed  discontinuity  at  an  abrupt 
boundary  can  be  significant.  Thus,  to  the  extent  that  the  bubble  clouds  have  well- 
defined  boundaries,  they  can  be  modeled  by  homogeneous  irregularities  with  an  internal 
sound  speed  given  by  (6).  One  finds  that  the  200-Hz  cross  section  of  a  meter-sized  object 
can  be  ss  20  dB  higher  than  the  corresponding  surface  backscatter  level  at  low  grazing 
angles.  For  example,  a  1-m  cylindrical  void  has  a  200-Hz  backscatter  cross  section  of 
~  —30  dB,  whereas  the  surface  backscatter  level  in  consistent  units  (see  Section  II)  is 
less  than  —40  dB  at  grazing  angles  below  20°.  Allowing  for  finite  permeability  does  not 
change  this  picture  significantly. 

Thus,  if  one  accepts  a  homogeneous  void  or  a  weakly  permeable  object  as  a  rea¬ 
sonable  first  approximation  to  a  subsurface  bubble  cloud,  the  problem  becomes  one  of 
determining  the  scattering  characteristics  of  such  an  object  near  an  irregular  surface. 
Although  highly  idealized,  this  model  can  identify  the  most  important  characteristics  of 
scattering  by  near- surface  scattering  objects.  It  is  possible  to  treat  this  composite  scat¬ 
tering  problem  by  generating  multiple  surface  boundary  integrals  that  enclose  each  ob¬ 
ject;  however,  this  approach  is  inefficient.  For  example,  a  method-of-moments  (MOM) 
computation  of  the  mutual  interaction  between  two  objects  of  comparable  size  requires 
an  eight-fold  increase  in  computation  time  over  that  required  for  a  single  object.  At 
Vista  Research,  Inc.,  we  have  developed  a  novel  method  that  allow-  us  to  use  the  known 
scattering  characteristics  of  the  individual  scatterers  in  an  exact  calculation  of  the  mu¬ 
tual  interaction  (multiple  scattering)  between  the  two  objects.  The  mutual-interaction 
computation  requires  less  time  than  the  computation  of  the  scattering  characteristics 
of  the  rough  surface  itself. 

The  method  is  completely  general  in  that  it  can  be  used  to  calculate  the  mutual 
interaction  between  any  pair  of  scattering  objects.  To  validate  the  method,  we  have 
performed  scatter  calculations  for  pairs  of  cylinders,  a  cylinder  near  a  smooth  surface, 
and  a  cylinder  near  a  rough  surface  (Section  III, 3).  In  all  cases,  we  observe  backscatter 
enhancements  in  excess  of  12  dB  over  the  isolated  backscatter  level  of  the  single  cylinder 
as  the  cylinder  is  moved  away  from  the  second  object.  The  coherent  superposition  of  the 
backscatter  from  two  identical  objects  or  an  object  and  its  image  can  account  for  a  6-dB 
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enhancement,  but  there  is  an  additional  pair  of  reciprocal  backscatter  paths  involving 
multiple  scatter  between  the  two  objects.  The  enhancement  is  somewhat  larger  for  a 
cylinder  adjacent  to  a  planar  or  rough  surface.  When  the  cylinder  is  in  contact  with 
a  plane,  the  scatter  cross  section  is  reduced  well  below  the  isolated  backscatter  level 
of  the  cylinder,  as  is  expected  from  image  theory.  Thus,  as  a  particle  breaks  away 
from  a  surface,  there  is  a  very  steep  rise  (%  30  dB  over  the  rough  surface  level)  in  the 
backscatter  cross  section. 

To  the  extent  that  this  mechanism  is  active  for  bubble  clouds  generated  by  breaking 
or  spilling  waves,  it  clearly  supports  the  bubble-cloud  hypothesis  as  an  explanation  of 
the  observed  anomalous  acoustic  backscatter.  The  large  backscatter  enhancement  asso¬ 
ciated  with  a  scattering  object  near  a  reflecting  surface  is  not  in  itself  a  new  result,  but 
there  are  many  ramifications  of  this  phenomenon  that  have  not  been  thoroughly  inves¬ 
tigated.  The  analytical  complexity  of  the  problem  has  been  the  principal  barrier,  but 
the  new  method  we  have  developed  reduces  the  computational  burden  to  a  manageable 
level. 


II  Background 

A  time-honored  problem  in  scattering  theory  is  the  computation  of  the  mutual  inter¬ 
action  among  disjoint  scattering  bodies.  Examples  of  this  type  of  problem  that  are 
important  in  active  acoustics  include  the  interaction  of  a  large  transducer  array  on  the 
ocean  bottom  and  the  scatter  from  near-surface  bubble  clouds.  The  problem  admits 
a  rigorous  formulation  in  terms  of  multiple  surface  integrals  over  the  scattering  ob¬ 
jects  that  can  be  evaluated  numerically.  This  approach,  however,  makes  no  use  of  the 
known  free-space  scattering  characteristics  of  the  individual  scatterers,  and  the  result¬ 
ing  numerical  methods  are  generally  inefficient  in  that  the  same  results  can  be  obtained 
with  less  computation  even  if  the  scattering  characteristics  of  the  individual  objects  are 
unknown  a  priori  and  must  be  computed  numerically. 

Because  considerable  effort  is  typically  expended  to  determine  the  scattering  char¬ 
acteristics  of  complex  bodies,  it  is  clearly  desirable  to  use  this  information  when  the 
scattering  characteristics  are  modified  by  the  presence  of  terrain,  water,  or  other  scat¬ 
tering  objects.  Traditionally,  this  multiple-scatter  problem  has  been  treated  by  diagram 
methods,  which  can  give  good  approximations  to  the  average  field  moments  for  random 
distributions  of  large  numbers  of  scatterers,  but  the  basic  equations  are  very  difficult  to 
evaluate  exactly.  As  a  simple  illustration  of  the  principle  involved,  recall  that  standing 
acoustic  waves  in  a  pipe  with  discontinuities  obey  the  transmission-line  equations  as  do 
waves  that  obey  the  more  general  equations  of  EM  theory.  One  can  solve  the  standing 
wave  problem  by  summing  the  multiplicity  of  reflected  waves,  but  the  transmission-line 
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equations  are  more  efficiently  solved  directly. 

It  is  not  widely  recognized  that  essentially  the  same  method  can  be  applied  to  waves 
scattered  by  multiple  objects.  The  modes  on  the  transmission  line  are  replaced  by 
spatial  Fourier  modes  propagating  forward  and  backward  relative  to  imaginary  planes 
separating  the  objects,  which  are  themselves  replaced  by  scattering  planes.  A  pair  of 
coupled  first-order  difference  equations  that  characterize  the  interaction  of  the  fields 
crossing  the  planes  can  be  solved  for  the  multiply-scattered  wave  fields  crossing  the 
reference  planes  in  terms  of  the  wave-wave  scattering  functions  that  formally  replace 
the  objects.  Once  these  fields  are  known,  the  backscattered  and/or  transmitted  waves 
are  easily  computed.  For  a  single  object  pair,  the  solution  is  particularly  simple.  In 
effect,  an  equation  of  about  the  same  complexity  as  the  equations  that  characterize 
the  individual  scatterers  must  be  solved.  Because  the  computation  time  varies  with 
the  cube  of  the  number  of  unknowns,  the  direct  approach  will  require  (21V)3,  whereas 
solving  the  mutual-interaction  equation  requires  only  N 3  additional  computations  once 
the  two  scattering  functions  are  known.  For  m  independent  objects,  the  computation 
time  varies  linearly  with  m  versus  (mlV)3  for  the  direct  method. 

The  mutual  interaction  equations  are  exact  to  the  extent  that  we  have  an  exact 
expression  for  the  wave- wave  scattering  function  for  each  object.  For  an  incident  plane 
wave  with  wave  normal  vector  k±(K'),  the  dyadic  scattering  function  Hn(k±(K),  k±(K')) 
characterizes  the  spectrum  of  scattered  plane- wave  components  indexed  by  their  wave 
normal  k±(K).  The  notation  emphasizes  the  fact  that  each  wave  vector  is  constrained 
to  have  amplitude  k  —  where  A  is  the  wavelength  in  the  surrounding  medium. 

Thus,  the  z  component  of  the  wave  normal  vector  is  related  to  the  magnitude  of  the 
transverse  component  K  by  the  relation 

kz{K)  =  ±Vk2  -  K2.  (7) 

When  K  >  k,  (7)  is  purely  imaginary,  and  the  corresponding  wave  is  said  to  be  evanes¬ 
cent.  These  evanescent  waves  must  be  properly  accommodated  in  the  dyadic  scattering 
function.  The  method  of  defining  the  scattering  functions  over  all  wave  numbers  is 
developed  in  the  appendix.  The  numerical  simulations  themselves,  which  are  described 
in  detail  in  Section  III,  yield  the  two-dimensional  Fourier  spectrum  of  the  total  wave 
field  near  the  scattering  surfaces.  The  scattered  wave  field  is  the  difference  between  the 
total  wave  field  and  the  wave  field  that  would  exist  in  the  absence  of  scatterers.  The 
far-zone  limit  of  any  propagating  wave  field  is  simply  related  to  the  two-dimensional 
Fourier  transform  of  the  corresponding  wave  field  evaluated  on  a  plane  near  the  scatter¬ 
ing  object.  Thus,  all  quantities  of  importance  are  readily  obtained  from  the  numerical 
simulations. 

It  is  impractical  to  perform  numerical  simulations  with  incident  plane  waves  for 
large  objects.  To  minimize  truncation  errors  from  edge  discontinuities,  it  is  necessary 
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to  use  a  narrow  beam.  Let  Ej(r)  represent  the  field  of  the  beam.  In  our  calculations 
we  normalize  the  beam  so  that  the  total  integrated  far-zone  intensity  is  unity.  The 
normalization  factor  has  the  units  of  length  or  area.  We  define  the  unitless  differential 
scattering  coefficients 


7(k 


«>ki)  =  I 


r|E(r)|2/(HE,|2) 

r2|E(r)|2/(A|E;|2) 


two-dimensional 

three-dimensional. 


(8) 


Note  that  7(k,,ki)  is  defined  in  terms  of  the  total  wave  field  and  integrates  over  all 
solid  angles  to  unity  for  a  lossless  scatterer.  These  definitions  are  commonly  used 
for  surface  scatter,  although  the  surface  cross-sectional  area  is  often  used  in  place  of 
the  cross-sectional  area  normal  to  the  beam.  In  [1],  we  distinguish  the  differential 
scattering  coefficient  7  from  the  differential  scattering  cross  section  by  using  the  notation 
a  =  7sin0i,  where  0;  is  the  incidence  angle.  For  particles,  however,  an  unnormalized 
cross  section  with  units  of  area  or  length  is  defined  in  terms  of  the  scattered  field  for  an 
incident  plane  wave.  The  standard  particle  cross  section  integrates  to  the  total  scattered 
intensity. 

The  highly  dissimilar  cross-section  measures  create  a  dilemma  for  data  presentation. 
We  have  chosen  to  resolve  it  by  retaining  the  definition  (8)  for  particles  as  well  as  sur¬ 
faces.  In  effect,  we  illuminate  the  particle  with  the  same  finite  beam  that  necessarily 
illuminates  the  surface.  When  presented  this  way,  the  particle  backscatter  coefficient  is 
in  units  that  can  be  compared  directly  to  the  surface  scatter  coefficient  or  the  cross  sec¬ 
tion  per  unit  area  if  the  scattering  coefficient  is  multiplied  by  the  cosine  of  the  incidence 
angle.  We  develop  formulas  to  convert  to  particle  cross-section  units  in  Section  III. 2, 
but  it  should  be  kept  in  mind  that  the  ideal  plane- wave  cross  section  is  a  limiting  form 
that  can,  in  a  numerical  simulation  or  a  real  measurement,  only  be  approached. 

Finally,  we  have  developed  our  analysis  in  its  most  general  form  for  vector  EM  fields. 
The  scalar  wave  fields  appropriate  to  the  acoustic  scatter  problem  are  easily  extracted. 
Indeed,  for  pressure  release  surfaces,  the  horizontally  polarized  EM  problem  and  the 
acoustic  problem  are  mathematically  identical,  as  we  showed  in  Appendix  A  of  [lj. 
One  need  only  interpret  the  square  root  of  the  relative  permittivity  as  a  velocity  ratio. 
The  parameters  chosen  for  all  our  examples  were  motivated  by  low-frequency  active 
acoustics.  We  take  the  velocity  of  sound  in  sea  water  as  1500  m/s. 
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Ill  The  Mutual  Interaction  Method 


III.l  Theory 

In  a  recent  paper  [10],  we  showed  that  the  mutual  interaction  among  many  scattering 
objects  can  be  characterized  exactly  by  a  pair  of  vector  wave  transport  equations.  This 
formulism  assumes  that  the  projections  of  these  scatters  onto  the  z  axis  do  not  overlap, 
whereby  the  medium  can  be  divided  into  slabs  bounded  by  planes  at  z  =  zn  such 
that  the  nth  object  lies  entirely  within  the  slab  zn_!  <  z  <  zn.  Each  discrete  object  is 
characterized  by  a  dyadic  scattering  functions  H^k*,  k'*),  where  the  first  wave  normal 
vector  refers  to  the  scattered  wave,  the  second  wave  normal  vector  to  the  incident  wave, 
and  the  sign  to  the  z  direction  of  propagation  for  real  kz  or  attenuation  for  imaginary 
kz.  Each  object  interacts  only  with  the  wave  fields  that  enter  the  slab  containing  it.  The 
exact  difference  equations  for  the  total  wave  fields  are  given  by  the  difference  equations 


fit  =  fit,  «?{»*»(*)£.}  +  Hr@Et,  +h;-®E; 


(9) 


and 


=  fi;  exp{sMK)£»}  +i£T@E;  +  K+ofij.,,  (10) 

where  K  =  |K|  is  the  transverse  wave  number,  Ln  =  \zn  —  zn_i\  is  the  width  of  the  slab 
containing  the  nth  object, 


and 


The  dyadic  scattering  functions  that  appear  in  (9)  and  (10)  are  the  natural  product 
of  integral  equation  computation  methods.  Formally,  they  represent  the  spectra  of  plane 
waves  scattered  by  each  incident  vector  plane  wave.  The  operator  @  is  defined  as 


E*  =  E±(K;  zr 

(ii) 

[1  -  K2lk2]1'2 

for  K  <  k, 

(12) 

i[K2/k 2  -  l]1/2 

otherwise. 

=  exp  {ikg(K)rj  jj  ^ 


dK' 
(2*)2 

ifcexp{— i(K  —  K')*pn} 


«U(K)-H°  (ku(K),  k±(K'))  (13) 


2  g(K) 

•En-i(K')  exp{ikg(K')i*}. 

n 

The  notation  E^_,  stands  for  E^_!  or  E~,  the  superscript  u  stands  for  +  or  £*  is  the 
distance  from  the  object  to  the  right  (  +  )  or  the  left  (  — )  bounding  plane  (see  Figure  2). 
pn  is  the  tranverse  position  vector  of  the  object,  and  k.u(K)  =  I  -  &-2kuku  where  I 
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is  a  unit  dyadic.  The  dyadic  £  is  a  projection  function  that  accounts  for  the  spatial 
dispersion  of  the  scattered  waves.  The  dyadic  is  the  scattering  function  of  the  nth 
object  when  it  is  located  at  the  origin.  The  factor  exp{— i(K  —  K')-pn}  accounts  for 
the  transverse  translation  from  the  origin. 

The  terms  on  the  right-hand  sides  of  (9)  and  (10)  admit  physical  interpretations. 
For  example,  there  are  three  contributions  to  the  total  forward-traveling  wave  at  the 
plane  z  =  zn.  From  (9),  we  see  that  the  first  term  represents  the  direct  propagation  of 
the  forward  wave  from  the  plane  z  =  zn_i  to  z  =  zn  as  though  the  medium  were  ho¬ 
mogeneous.  The  second  term  represents  the  incremental  contribution  due  to  scattering 
of  the  forward  wave  in  the  forward  direction.  The  last  term  represents  the  incremental 
contribution  from  the  backward  wave,  which  propagates  from  the  plane  z  —  zn  to  the 
object  and  scatters  in  the  forward  direction.  A  similar  interpretation  applies  to  the 
corresponding  terms  in  (10). 

Equation  (13)  also  admits  a  physical  interpretation.  The  last  two  factors  in  the 
integrand  represent  the  incoming  wave  as  it  would  be  seen  at  the  plane  passing  through 
the  center  of  the  object  if  it  were  to  propagate  freely  from  the  appropriate  boundary 
plane.  The  object,  now  acung  only  on  the  field  at  the  center  plane,  scatters  the  incoming 
wave,  which  produces  the  incremental  scattered  waves  represented  by  the  dot  product 
between  the  incoming  wave  and  the  scattering  function  weighted  by  ik/2g(K).  These 
scattered  waves  then  propagate  freely  to  the  slab  boundary.  The  integral  operation 
accounts  for  the  cumulative  scattering  from  all  Fourier  components  of  the  incoming 
wave.  Thus,  in  the  spectral-domain  formulism  the  scattering  objects  are  effectively 
replaced  by  scattering  planes.  The  appendix  discusses  in  detail  how  to  calculate  the 
scattering  function  in  a  coordinate  system  centered  on  the  object  from  the  spectrum  of 
waves  scattered  by  an  incident  vector  plane  wave. 

The  utility  of  equations  (9)  and  (10)  for  solving  complex  problems  can  be  illustrated 
by  considering  their  application  to  the  computation  of  the  exact  mutual  interaction 
between  two  arbitrary  scattering  objects.  Solving  the  two-object  problem  using  MOM 
requires  an  eight-fold  increase  in  computation  time  and  twice  the  storage  as  is  needed 
for  a  single  object.  Twersky  [11]  has  formulated  the  two-object  multiple-scatter  problem 
as  an  infinite  multiple-scatter  series,  each  term  of  which  uses  the  asymptotic  scattering 
dyadic  for  one  or  more  of  the  particles.  His  method  yields  some  useful  approximations 
for  random  particle  configurations,  but  the  multiple-scatter  series  cannot  be  solved 
exactly.  On  the  other  hand,  by  simply  writing  (9)  and  (10)  for  a  two-particle  system, 
we  obtain  an  exact  linear  equation  for  the  fields  crossing  the  plane  separating  the  two 
scatterers.  Once  these  excitation  fields  are  computed,  the  complete  problem  can  be 
solved  with  a  few  simple  manipulations;  moreover,  one  scatterer  can  be  replaced  with 
a  reflecting  plane  or  a  rough  surface. 

For  two  objects  separated  by  a  distance  AL  as  shown  in  Figure  2,  a  straightforward 
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Figure  2:  Geometry  for  two-object  computation  with  definition  of  symbols. 


application  of  (9)  and  (10)  leads  to  the  following  equation  for  Ej  : 

Er -S2-+@(lr@Er)  =  %nc**p{ikgL3}  +fi2-OEine>  (14) 

where  Ejnc  is  evaluated  at  the  plane  z  —  z2.  This  equation  can  be  solved  numerically 
solved  for  E^.  The  remaining  fields  are  then  determined  from  the  following  equations: 

E+  =  Sr@Er  (15) 

E„  =  Ef  exp{zfc^Li} +Iii_~@Ef  (16) 

E+  =  E+exp{%L2}+S2++@E++E2+-@Einc  (17) 

These  results  are  valid  for  any  pair  of  objects,  although  some  care  must  be  taken  in 

defining  the  appropriate  scattering  function  for  a  conducting  surface  of  infinite  extent. 

Suppose,  for  example,  that  the  first  object  is  a  perfectly  conducting  infinite  planar 
surface.  An  arbitrary  surface  can  be  generated  by  adding  higher-order  spatial  Fourier 
components  to  the  zeroth-order  term  representing  the  mean  surface  height.  This  com¬ 
ponent  scatters  the  incident  wave  in  the  specular  direction  and  cancels  the  incident  wave 
in  the  forward  direction  on  the  other  side  of  the  surface.  Thus,  the  scattering  function 
for  an  infinite,  perfectly  conducting  surface  contains  a  delta  function.  The  scattering 
function  for  a  conducting  plane  is  given  by  (.475,  .476)  in  the  appendix.  It  is  convenient 
to  decompose  the  scattering  function  into  a  regular  part.  and  a  singular  part.  The 
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singular  part  consists  only  of  a  delta  function,  which  can  be  integrated  analytically. 
The  four  scattering  functions  for  a  perfectly  conducting  surface  can  be  written  formally 
as 

fiHK.K')  =  j^-I^K.K')  - -  K'XX  -  2a*az),  (18) 

and 

*±-fi1±:t(KlK')  =  iS±-fif:tt(KfK/)  -  -  K%  (19) 

where  a*  is  the  unit  vector  along  the  z  direction.  Upon  substituting  the  above  scattering 
functions  into  equations  (14)  through  (17),  we  obtain  the  following  results  for  a  surface- 
particle  aggregate: 

Er  -Ej^dr^Er)  -|fc+@[(2a,a, -iyE;(K')exp{i2kg(K')et}} 

=  Eincexp{ikgL2}  @Einc  (20) 


E+(K)  =  (2a,ai -I)*Er(K)exp{i2fc5(KK  +  } +  H+-t(K,K')@E-(K')  (21) 

Eo(K)  =  0  _  (22) 

E+(K)  =  Et(K)exp{ikg(K)L2}+llt+@Et +}&-@Einc  (23) 


To  convert  the  forward  and  backward  wave  fields  to  physical  observables,  we  note 
that  at  distances  well  removed  from  the  reference  plane,  the  Helmholtz  integral  can  be 
approximated  by  the  method  of  stationary  phase  to  yield 

{ikg(Kr)E:k(Kr)  exp{ikr}  /  (2irr)  three-dimensional 

(24) 

\J — ik  g{KT)E±{'KT)exp{ikr}  j  two-dimensional 


where  Kr  is  a  wave  vector  directed  along  r.  The  physical  observable  of  particular 
interest  to  us  is  the  angular  power  density  in  the  direction  k  for  a  given  incidence  k;. 
It  is  given  by 


P±(k(K),kl) 


K’2|gf(K)j2|E±(K)|2/(47r2)  three-dimensional 
,  fc|3(A’)|2|E:t(K)|2/(27r)  two-dimensional 


(25) 


In  the  next  section  we  shall  relate  the  angular  power  density  to  the  dimensionless 
differential  cross  section  7(k(K),k;). 

It  is  instructive  to  specialize  the  above  results  for  the  case  of  parallel  polarization. 
In  this  case  the  wave  fields  become  scalar  and  equations  (14)  through  (17)  or  (2d) 
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through  (23)  can  be  converted  to  their  equivalent  scalar  forms  by  making  the  following 
substitutions: 


fij  -  a 

«”(K).H:(lc“(K),k”(K')).£:.,(K')  -  ir;(k“(K),k'(K'))*-i(K') 
(2a,a,-I).£r(K)  . - fr(K) 


The  equations  are  greatly  simplified  for  isotropic  scatterers,  in  which  case 

*“(K).|£(ku(K), ku(K'))  -  An,  (26) 

where  hn  is  a  constant.  When  (26)  holds,  the  @  operator  can  be  evaluated  numerically. 
For  example,  (13)  becomes 

Br®£:-  -  r4^7rknexp{— t[K-pn  -  kg{K)l„  ]}^-i(pn.  V>).  (2?) 

n  ^9{K) 

where  z™  is  the  z  coordinate  of  the  nth  particle  and  ipn-i  is  the  free-space  propagation 

§  n 

of  the  wave  field  from  the  slab  boundary  zn_i  or  zn,  i.e., 


*£-(/>„>  *.)  =  JJ  T^ps  ^-.(K)eap{tiK.a„  +  kj(K)Up,-2.-.|]}  (28) 


(2tt) 


With  these  simplifications,  the  spatial-domain  forms  of  (14)  or  (20)  can  be  computed 
analytically.  To  accommodate  both  the  two-  and  three-dimensional  forms,  we  let 


rr  dK  iexp{i[K-p  +  kg(K)\z\}}_ 

=  JJ  (2 2 kg(K) 

exp{iAr}/(47rr)  three-dimensional 

=  -  (29) 

i#o^(fcr)/4  two-dimensional 


Note  that  the  three-dimensional  form  has  the  units  m  1 ,  whereas  the  two-dimensional 
form  has  no  units. 

3y  using  (27)  and  (29),  the  scalar  isotropic  form  of  (14)  can  be  solved  for  l,  zp i)- 
To  this  end,  we  multiply  both  sides  of  (14)  by  exp {ikg{K)£^},  take  the  inverse  Fourier 
transform,  and  evaluate  the  result  at  p  =  pl  to  obtain 


i  (pi.zpi) 


(Pi,  zpi )  +  h2k2G  12  ^Pinc  (p  i.  fpj) 
1  -  h1hi(k2Gl2)2 


(30) 


where 

G 12  =  G{\pi  —  f>2 


(31) 


« 


< 
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is  the  spherical-wave  propagation  factor  for  the  distance  between  the  particles.  The 
denominator  of  (30)  accounts  for  the  mutual  interaction  of  the  particles.  From  (17), 
the  backscattered  wave  spectrum  can  be  computed  as 


tfc(K)  =  ifc[20(/C)]-1/iiexp{-i[K-pi  -  kg(K)(AL  +  1$)}}  (px,  zpl) 

+  ik3[2g(K))~lhih2Gi2  exp{— »[K*p2  -  kg(K)tf]}  ^ {px,  zpl) 

#  +  ik[2g(K)]~1h2exp{~i[K‘p2  -  kg(K)£2}}tpinc(p2,zp2).  (32) 

The  three  terms  in  (32)  represent  the  direct  propagation  of  %l)x ,  the  forward  scattering 
of  ifii  by  the  second  object,  and  finally  the  incident  wave  backscattered  by  the  second 
object. 

^  To  compute  the  forward-scattered  wave  field  spectrum  i/>q(K),  however,  we  must 

proceed  in  two  steps.  From  (14), 

dx(K)  =  xf;inc(K)exp{ikg(K)L2}  +  ik[2g{K)}~lh2  exp{-i[K-p2  -  kg{K)l2}} 

*  [Vw(p2,  zp2)  +  h1k2G12ip;(p1,zpl)Y  (33) 

®  Now  t/;0-(K)  can  be  computed  from  (16)  as 

4>o(K)  =  $i{K)exp{ikg(K)Lx}  +  ik[2g{K)}~1  hx 

xexp{-i(K*p1  -  kg(K)£x\}  ipx(px,zpX).  (34) 


Following  the  same  procedure,  equation  (20)  can  be  solved  for  ipx{p2,  -A L)  for  the 
case  of  an  isotropic  cylinder  above  a  perfectly  conducting  plane.  However,  the  factor 
used  is  exp{ikg( K)(/\L  +  £  f )}  and  the  inverse  Fourier  transform  is  evaluated  at  p  =  p2. 
Explicitly, 


A  (p2,-AL) 


^mc(p2,  -  A L)  +  h7k2  G(2A£)  U>inc(p2,  -Zpz) 
1  4-  h2k2  G(2AL) 


The  backscattered  wave  field  is  given  by 


(35) 


&(K)  =  —  V^i"(K) exp{ikg(K)(2i  x  4-  L2)}  +  ik{2g(K)}~x h2 

xexp{-i[K-p2  -  kg(K)£2  ]}  \il>inc{p2,  zp2)  -  t/q- (p2,  -AL)] ,  (36) 

where  xpx  (K)  is  found  from  (20)  to  be 


fr(K) 


4>inc{K)exp{ikg(K)L2}  +  ik\2g{K)\  lh2 

x  exp{-i[K-p2  -  kg(  K)f  7  }  rtnc(p2.  cp2 )  -  c'1"(p2. -AL)| .  (37) 
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We  shall  see  that  the  two  terms  in  (30)  or  (35)  combined  with  the  leading  two 
terms  in  (32)  or  (36)  can  add  coherently  to  produce  a  four-fold  amplitude  enhance¬ 
ment  or  12  dB.  For  irregular  surfaces,  which  cannot  be  accommodated  analytically,  the 
enhancement  is  potentially  larger. 

In  conclusion  let  us  note  that  if  the  scattering  functions  are  known,  (14)  or  (20)  can 
be  put  in  the  form  of  a  system  of  linear  equations  by  performing  a  matrix  multiply  to 
evaluate  the  left-hand  sides  and  a  vector  multiply  to  evaluate  the  right-hand  sides.  The 
system  of  equations  can  then  be  solved  numerically  to  determine  Ej\  The  computation 
time  is  dominated  by  the  matrix  multiply  (~  N?  computations  if  H(k*,kf)  is  an 
Nr  x  Nr  matrix),  and  the  solution  of  the  resulting  matrix  equation  (an  additional 
~  N?  computations).  The  computation  time  for  the  remaining  fields,  which  involve 
only  matrix-vector  multiplies,  is  negligible.  Thus,  the  total  number  of  computations  for 
solving  the  two-scatterer  problem  varies  like  2 iVT3  with  the  mutual-interaction-method 
(MIM)  versus  ( 2Nr )3  for  MOM. 
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III. 2  Numerical  Simulations 


For  the  numerical  simulations,  we  will  only  consider  problems  that  admit  variations  in 
two  dimensions  ( x,z ).  The  above  equations  are  essentially  unchanged  except  that  K  is 
replaced  by  kx.  For  example,  equation  (13)  becomes 


H^OEt,  £  exp {ikg(k,)rj  J 


dk'x  ik  exp{— i(kx  —  A£)xn} 


2 g{kx) 


To  solve  for  the  mutual  interaction  between  two  objects,  we  need  to  know  their  individ¬ 
ual  scattering  functions.  The  appendix  presents  a  complete  derivation  of  the  scattering 
function  of  a  single  object  with  explicit  results  for  an  infinite  cylinder,  an  infinite  con¬ 
ducting  plane,  and  a  planar  conducting  surface.  The  scattering  function  of  all  but  the 
rough  conducting  surface  can  be  obtained  in  analytic  form. 

To  illustrate  the  mutual  interaction  method,  we  will  perform  computations  for  (1) 
a  pair  of  parallel  cylinders  and  (2)  a  cylinder  above  an  infinite  conducting  surface.  The 
z  axis  is  vertical,  with  the  axis  of  cylindrical  symmetry  lying  along  the  y  axis.  Thus, 
the  surface  height  varies  only  in  the  x  direction.  The  incident  wave  normal  vector  is 
assumed  to  lie  entirely  in  the  xz  plane  with  two  principal  polarizations.  For  axial  or 
parallel  polarization,  the  incident  electric  field  vector  is  perpendicular  to  the  incidence 
plane  (i.e.,  parallel  to  the  cylinder  axis).  For  transverse  or  vertical  polarization  it  lies 
in  the  plane  of  incidence  (i.e.,  perpendicular  to  the  y  axis).  The  projection  of  the 
scattering  dyadic  onto  the  direction  of  the  incident  wave  can  be  written  in  the  general 
form 

^•IJ( A:JX,  A;,x),at( A),x)  =  ^s{ksx)H\\j^{ksx,  klx),  (39) 

where  al<s(kx)  is  a  unit  vector  along  the  kx  Fourier  component  of  the  incident  or  scattered 
wave  field.  In  general,  a,  depends  on  the  object’s  geometry  and  the  incident  wave 
polarization. 

For  an  infinite  cylinder  with  refractive  index  m  and  radius  a,  lying  along  the  y  axis, 
the  scattering  functions  are  given  as  follows: 


OO 

#||(k,,  k.)  =  Mk~2  cn  exp(ma),  (40) 

71=  —  OO 


and 

OO 

tf±(k,,  k.)  =  4ik~2  dn  exp(ina),  (41) 

71=  -  OC 

where  cn,  dn,  and  a  are  given  by  (,463)  to  (.472)  in  the  appendix.  These  equations 
define  the  scattering  functions  for  all  values  of  ksx  and  klx.  When  \ksx\  k  or  \klxl  «  k. 
a  is  complex. 
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Figure  3  shows  the  magnitude  of  the  scattering  functions  of  a  perfectly  conducting 
cylinder  and  a  dielectric  cylinder  with  refractive  index  m  =  1.5  illuminated  by  an 
axially  polarized,  7.5-m  free-space  wavelength  wave.  Both  cylinders  have  the  same  1-m 
diameter.  The  plot  shows  the  magnitude  as  a  function  of  the  normalized  scattered 
wave  number  k,x  for  a  fixed  incidence  k{x  =  0.  The  scattering  functions  are  essentially 
constant  over  the  range  of  \ktx/k\  <  1,  i.e.,  all  real  azimuthal  angles.  This  is  the 
scattering  behavior  of  the  zeroth-order  cylindrical  normal  mode,  which  is  azimuthally 
independent.  Thus,  the  zeroth-order  cylindrical  normal  mode  dominates  the  scattering 
for  this  cylinder. 

In  general,  a  small  cylinder  (a/ X  <  1)  scatters  essentially  uniformly  in  azimuth,  and 
the  scattering  function  can  be  approximated  by  the  zeroth-order  term  only,  whereby 

{  h[  }  =  4tk  2  {  "  h’  ^ 

which  can  be  used  in  (30)  and  (32)  or  (35)  and  (36).  If  the  cylinder  diameter  is  large 
compared  to  the  incidence  wavelength,  however,  the  scattering  function  is  no  longer 
axially  symmetric.  This  behavior  can  be  seen  in  Figure  4,  which  is  the  corresponding 
scattering  function  for  a  large  cylinder.  In  this  case  the  higher-order  cylindrical  normal 
modes  are  significant,  and  the  integral  equation  (14)  lends  itself  to  numerical  solution 
only.  For  rough  surfaces,  numerical  solutions  must  be  used  in  all  cases. 

The  scattering  functions  for  an  infinite  conducting  plane  are  given  by  (>175)  and 
(i476),  which  are  reproduced  here  for  reference  purposes: 

kix)  =  -  kx)(ax a*  -  aza2)  (43) 

k±-H±±(^,*,^,i)  =  - s(k>*  ~  fci*)(a*a*  +  a*az)  (44) 

Thus,  the  scattering  function  for  a  conducting jdane  consists  only  of  a  delta  function. 
The  backward  scattering  functions  H+~  and  H~+  describe  the  specular  reflection  of 
the  incident  wave.  The  forward  scattering  functions  H++  and  H--  represent  the  di¬ 
rect  propagation  of  the  incident  wave  with  a  phase  shift  of  180°  so  that  the  total  wave 
vanishes  identically  behind  the  plane.  The  scattered  wave  behind  a  flat  or  rough  con¬ 
ducting  surface  always  exactly  cancels  the  incident  wave.  Hence,  the  forward  scattering 
function  always  contains  a  delta  function.  The  scattering  function  can  have  additional 
delta  functions,  however,  because  the  scattering  of  a  wave  from  a  surface  is  a  wave- wave 
scattering  process  where  the  second  wave  comes  from  the  surface  itself.  The  spectrum 
of  a  rough  surface  can  contain  one  or  more  discrete  Fourier  components.  The  discrete 
component  at  zero  spatial  frequency  corresponds  to  the  zero- mean  surface  component. 
These  discrete  Fourier  components  are  responsible  for  the  specular  reflection  of  the 
incident  wave. 
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Scattering  Function  H  +  +  (Axial  Polarization) 


Normalized  Scattered  Wave  Number  ksx/k 


Scattering  Function  FH —  (Axial  Polarization) 


Normalized  Scattered  Wave  Number  ksx/kO 


Figure  3:  Magnitude  of  the  scattering  function  of  an  infinitely  long  cylinder  15  m  in 
diameter  illuminated  by  an  axially  polarized  7.5-m  wavelength  wave  incident  along  the 
z  axis  for  two  cases:  (a)  perfect  conductor.  (6)  dielectric  with  refractive  index  of  1.5. 


Scattering  Function  H  +  +  (Axial  Polarization) 


Scattering  Function  H  +  -  (Axial  Polarization) 


Figure  4:  Magnitude  of  the  scattering  function  of  an  infinitely  long  cylinder  15  m  in 
diameter  illuminated  by  an  axially  polarized  7.5-m  wavelength  wave  incident  along  the 
z  axis  for  two  cases:  (a)  perfect  conductor,  (6)  dielectric  with  refractive  index  of  1.5. 
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Because  each  specular  reflection  is  characterized  by  a  delta  function,  it  is  desirable  to 
separate  the  scattering  function  into  a  regular  part,  H*,  and  a  singular  part  consisting 
of  all  the  delta  functions.  Thus,  the  scattering  dyadic  can  be  written  formally  as 

H^(^«*)^tx)  "b  h'pqS(kax  kXfpq ) ,  (4b) 

PA 

where  hpq  is  the  appropriate  amplitude,  ap  or  a,  is  a  unit  vector,  and  the  summation 

runs  over  all  three  orthogonal  directions.  The  singular  part  can  then  be  integrated  ana¬ 

lytically  in  the  above  difference  equations  before  any  numerical  solutions  are  attempted. 

The  above  decomposition  is  always  applied  to  the  surface  scattering  functions  in  ov r 
numerical  solutions.  For  a  conducting  surface  the  regular  part  is  identically  zero.  It  is 
not  obvious,  however,  how  to  identify  the  singular  part  for  a  rough  surface  since  it  is 
implicitly  contained  in  the  numerical  solutions.  In  fact,  the  results  from  (480)  to  (A95) 
in  the  appendix  can  be  summarized  as  follows: 

H^~(k.xikix)  =  Mh^lE.y(k.x)  (46) 

H+~(k.x,kix)  =  )  (47) 

H{J-(k.x,kix)  =  HZ~(k.x,kix)  =  --9£3x) 6(k,x  -  klx)  (48) 

To  extract  the  singular  part,  note  that  the  backward-scattering  function  H+~  is  always 
singular  because  the  zero-mean  component  of  an  infinite  rough  surface  always  scatters 
the  incident  wave  in  the  specular  reflection  direction.  Assuming  that  the  surface  has  no 
other  discrete  harmonics,  we  can  identify  the  regular  part  of  the  backward  scattering 
functions  by  subtracting  the  flat-surface  component  as  follows: 

H,u'  =  -  *»)  («) 

The  above  results  are  given  in  terms  of  the  scattered  wave  fields  E ,  and  E,y,  but 
these  fields  cannot  be  derived  analytically  for  a  rough  surface.  In  principle,  they  can 
be  obtained  by  solving  numerically  the  problem  of  a  unit- amplitude  incident  plane 
wave  scattering  from  the  surface.  Unfortunately,  the  incident  wave  spectrum  is  a  delta 
function  and  therefore  cannot  be  implemented  numerically;  however,  the  delta  function 
can  be  approximated  by  a  finite-amplitude  function  to  as  high  a  degree  of  accuracy  as 
desired.  A  convenient  choice  is  a  narrow  gaussian  incident  beam  of  the  form 

E,{kx]  z)  =  exp{-(/cr  -  kxl)2 /kl}exp{ik:(kxl)z},  (50) 
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where  k =  ksindi  and  kw  =  const.  (This  gaussian  beam  is  also  used  in  the  appendix 
in  the  discussion  of  analytic  continuation  of  the  scattering  function.)  Note  that  the 
fc*-integrated  magnitude  is  unity,  i.e., 

/I  lr  =  ( 5i) 

and  the  incident  beam  approaches  a  unit-amplitude  plane  wave,  27T 6(kx  —  kxi),  as  kw  — * 
0.  The  scattering  function  is  very  sensitive  to  the  beamwidth  parameter  kw.  This 
parameter  is  chosen  such  that  the  numerically  discretized  gaussian  beam  has  only  one 
overwhelmingly  dominant  component  at  the  central  frequency  kx  =  kx{.  We  typically 
choose  kw  =  8 /  L  where  L  is  the  length  of  the  discretized  surface. 

The  mutual  interaction  integral  equation  (14)  is  readily  solved  after  the  scattering 
function  of  each  scatterer  is  completely  determined,  either  numerically  or  analytically. 
The  scattered  wave  field  Ej“  can  be  solved  numerically  for  an  arbitrary  incident  wave 
field,  from  which  Ej  and  all  the  backscatter  attributes  can  be  computed  (see  Figure 
2);  however,  these  backscatter  attributes  are  usually  defined  for  an  incident  plane  wave. 
For  numerical  computation,  we  can  approximate  a  plane  wave  by  a  narrow  gaussian 
beam  of  the  form 


Ei(kx;z)  =  ^exp{-(^  -  A«-)V*i}exp{*Mki)*}>  (52) 

where  A  is  the  amplitude  to  be  chosen.  The  choices  were  discussed  in  Section  II. 

For  sufficiently  small  values  of  kw/k ,  which  approximates  a  plane  wave,  it  can  be 
shown  that  if  .4  is  chosen  so  that  the  amplitude  of  the  narrow  gaussian  beam  integrates 
to  unity,  then  the  angular  power  density  P(k,,k;)  of  the  scattered  wave  field  in  the 
direction  k,  [cf.  Eq.  (25)]  is  numerically  equal  to  the  cross  section  o-(k,,kt)  of  the  two- 
object  aggregate.  If  .4  is  chosen  so  that  the  power  in  the  beam  integrates  to  unity,  then 
P(k4,ki)  is  equal  to  the  differential  scattering  coefficient  7(k,,k^).  In  the  former  case, 
the  scaling  factor  A  is  equal  to  2y/Tr/kw  (cf.  Eq.  (50)],  whereas  in  the  latter  case  .4  is 
determined  by  the  condition  that  the  total  power  integrated  over  all  directions  is  unity, 
ie., 

[lit  rk  rjU  ^ 

l  |Ei(r,0)|2r  d6  —  J  ^  g(kx)\Ei(kx]  r)|2  =  1.  (53) 

Thus,  the  conversion  factor  in  two-dimensional  problems  is  given  by 


g~(k,,ki) 

7(k.,ki) 


£7  J_k  ^0(**)exp{— (fc,  -  kxt)*/kl}. 


(54) 


For  higher-than-grazing  incidence  this  integral  can  be  evaluated  in  a  power  series  as 
follows: 


(7( k# ,  k, )  g( kxt )  ^  r  i  ~n  1  \ ;  2 
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Figure  5:  Backscatter  conversion  factor  cr0/70  for  A  =  7.5  m  as  a  function  of  incidence 
angle  for  different  values  of  the  normalized  beamwidth  parameter  kw/k. 


a0  =  1, 

(56) 

a2  =  —ll[2k2g2(kxi)\, 

(57) 

1  +  4  k2Jk2 

°4  8 fcV(*«)  ’ 

(58) 

1  +  I2k-2kit  +  8k-4ki- 

(59) 

_  Xt  1  XI 

1  Gk6g6(kxl) 

Note  that  cr  has  the  dimension  of  length,  whereas  7  is  dimensionless,  which  is  consistent 
with  (54).  Equation  (54)  was  calculated  numerically  for  A  =  7.5  m  and  is  plotted  in 
Figure  5  as  a  function  of  the  incidence  angle  for  different  values  of  the  normalized 
beamwidth  kw/k.  It  can  be  seen  that  for  the  same  beamwidth  ku..  the  conversion  factor 
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depends  on  the  angle  of  incidence.  All  of  the  examples  in  this  report  use  A  =  7.5  m 
and  kw/k  =  9.095  x  10-3  (which  corresponds  to  Lj\  =  140)  so  that  the  conversion 
factor  is  22.16  dB  at  60°  incidence,  20.51  dB  at  70°,  and  17.51  dB  at  80°.  Since 
crj 7  is  proportional  to  A2,  Figure  5  can  also  be  used  to  infer  the  conversion  factor  for 
wavelengths  other  than  7.5  m. 


III. 3  Examples 


To  illustrate  the  flexibility  of  MIM  for  different  objects,  the  calculation  procedure  de¬ 
scribed  in  Section  III.2  has  been  applied  to  the  following  scattering  problems:  (a)  two 
identical  closely  spaced  cylinders,  (b)  an  infinite  cylinder  near  a  conducting  plane,  and 
(c)  an  infinite  cylinder  near  a  conducting  rough  surface  with  a  power-law  spatial  wave 
number  distribution.  The  cylinders  arc  1  m  in  diameter  and  perfectly  conducting  (PC). 
The  wavelength  for  all  our  computations  is  7.5  m,  corresponding  to  an  acoustic  fre¬ 
quency  of  200  Hz.  We  use  the  scalar/ acoustic  model,  which  is  equivalent  to  parallel 
EM  polarization  for  the  two-dimensional  problem.  For  the  initial  calculations,  we  use 
an  incidence  angle  of  60°  with  respect  to  the  two-particle  vertical  axis  or  the  surface 
normal.  The  small  incidence  angle  is  the  least  demanding  from  a  computational  point 
of  view.  Thus,  it  is  a  convenient  starting  point. 


We  first  show  the  effects  of  the  strong  mutual  interaction  between  two  cylinders — 
example  (a).  Figures  6  shows  the  bistatic  scattering  coefficient  for  A L  —  5  m  super¬ 
imposed  on  the  scattering  cross  section  for  a  single  cylinder.  The  scattered  field  from 
the  single  cylinder  is  nearly  isotropic  at  a  level  —29.78  dB  in  the  normalized  cross  sec¬ 
tion  units  we  are  using.  The  60°  incidence  angle  with  respect  to  the  two-cylinder  axis 
corresponds  to  <j>  =  —45°  and  0fnc  =  15°  in  Figure  2.  Figure  7  shows  the  backscatter 
variation  as  the  cylinder  is  moved  away  from  the  surface.  The  backscatter  is  alternately 
enhanced  and  reduced  relative  to  the  free-space  level;  however,  the  peak  backscatter 
enhancements  exceed  the  6-dB  level  that  would  result  from  a  simple  coherent  super¬ 
position  of  the  individual  backscattered  fields.  Thus,  multiple  scattering  is  important 
here.  The  solid  curve  is  the  analytic  result  for  isotropic  scatterers — (30)  and  (32).  The 
dominant  features  are  preserved,  but  there  are  significant  differences  even  though  the 
departures  from  isotropy  are  slight.  We  shall  see  that  these  effects  are  more  pronounced 
for  plane  surfaces. 


To  validate  the  MIM  computations,  the  two-cylinder  MIM  results  have  been  com¬ 
pared  to  computations  made  with  the  NEC  generalized  electromagnetics  code.  The 
NEC  code  performs  three-dimensional  computations  using  MOM.  Figure  8  shows  the 
backscatter  level  as  a  function  of  separation  for  a  wave  incident  along  the  axis  of  the 
cylinders  ( <f>  =  O°;0inc  =  0°),  which  tends  to  miminimize  the  backscatter  enhancement. 
Figure  9  shows  the  same  computation  for  an  incident  wave  broadside  to  the  axis  <<f  the 
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%  Figure  6:  Differential  scatter  coefficient  for  (a)  a  single  cylinder  and  (6)  two  identical 

cylinders  separated  by  5  m. 


cylinders  ((f)  =  -45°;  0{nc  =  45°),  which  tends  to  maximize  the  backscatter  enhancement. 
As  discussed  in  Section  III.  1 ,  the  12-dB  enhancement  over  the  free-space  backscatter 
level  is  attributed  to  the  fact  that  the  small  cylinders  are  nearly  isotropic  radiators 
and  there  is  a  reciprocal  pair  of  paths  involving  multiple  reflections  that  can  contribute 
an  additional  doubling  of  the  backscatter  amplitude.  We  shall  see  that  this  effect  is 
even  more  prominent  for  small  cylinders  near  a  plane  or  rough  surface.  The  small 
discrepancies  between  the  MOM  and  MIM  results  are  attributable  to  the  three-  and 
two-dimensional  geometries  and  the  finite  beamwidth  used  in  the  mutual-interaction 
computations.  In  the  NEC  code,  an  incident  plane  wave  and  long  but  finite  cylinders 
were  used.  The  MIM  calculations  ran  more  than  three  times  faster  than  the  NEC  code, 
but  the  comparison  is  not  definitive  because  of  the  generality  of  the  NEC  code  and  the 
fact  that  the  MIM  computation  was  initiated  with  known  scattering  functions. 
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Two  PC  Cylinders  (o=-45  deg,  3^=15  deg) 


Figure  7:  Backscatter  cross  section  for  two  identical  cylinders  as  a  function  of  their 
separation  A  L. 
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Two  Perfectly  Conducting  Cylinoers  —  Bockscotter  -  90  deg 


Figure  9:  Comparison  of  NEC  and  MIM  calculations  for  a  wave  incident  broadside  to 
the  axis  of  the  cylinders. 
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20. 


Figure  10:  Differential  scatter  coefficient  for  a  cylinder  and  a  flat  surface  for  three  values 
of  the  distance  A L  (4,  5,  and  6  m). 

We  now  turn  to  example  (b),  effectively  by  replacing  the  first  object  in  Figure  2 
with  a  flat  surface  at  z  =  0.  Figure  10  shows  the  differential  scatter  coefficient  as  a 
function  of  scattering  angle  for  an  incidence  angle  of  60°  with  respect  to  the  surface 
normal  and  for  three  values  of  A L  (4,  5,  and  6  m).  The  peaks  at  6  —  -60°  are  due 
to  the  specular  reflection  from  the  surface.  In  the  backscatter  directions,  the  scatter 
is  much  stronger  than  the  corresponding  level  from  a  single  cylinder;  moreover,  the 
backscatter  at  A L  =  4  m  exceeds  the  maximum  level  achieved  for  two  cylinders — see 
Figure  7.  Figure  11  shows  the  variation  in  the  differential  scatter  at  A L  =  5  m  for 
incidence  angles  of  60°,  70°,  and  80°.  For  this  particular  separation,  the  backscatter 
level  is  nearly  constant  for  the  three  incidence  angles  shown,  but  the  incidence-angle 
dependence  of  the  backscatter  level  is  strongly  dependent  on  the  distance  of  the  cylinder 
from  the  surface. 
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Figure  11:  Differential  scatter  coefficient  for  a  cylinder  and  a  flat  surface  for  three  values 
of  the  incidence  angle  0,  (60°,  70°,  and  80°), 


Figures  12,  13,  and  14  show  the  backscatter  variation  as  a  function  of  distance  for 
the  60°,  70°,  and  80°  incidence  angles  used  in  Figure  11.  The  backscatter  enhancement 
peak  is  insensitive  to  incidence  angle,  but  its  location  moves  progressively  away  from 
the  surface  as  the  incidence  angle  is  increased.  For  comparison  purposes  we  have  also 
plotted  the  analytic  isotropic  formula  (36)  where  the  cylinder’s  scattering  function  is 
approximated  by  the  azimuthally  independent  term  only.  If  we  use  strictly  isotropic 
scattering  functions  in  the  numerical  computations,  the  analytic  and  numerical  results 
are  identical.  Thus,  the  differences  between  the  two  sets  of  computations  shown  in 
Figures  12,  13,  and  14,  are  due  entirely  to  the  small  departure  from  isotropy  of  the 
scattering  from  the  1-m  cylinder.  We  see  that  the  isotropic  formula  predicts  higher 
peaks  and  deeper  nulls,  with  the  differences  increasing  with  increasing  incidence  angle. 
Both  the  analytic  and  numerical  results  show  the  effects  of  higher-order  interaction 
terms  as  the  incidence  angle  increases.  Both  results  also  show  that  the  backscatter 
enhancement  peak  increases  with  increasing  incidence  angle.  At  80°,  it  exceeds  12  dB 
for  both  the  exact  and  isotropic  curves. 
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GammaO 


PC  Cylinder  Above  Conducting  Plane  (60  deg  Inc) 


Figure  12:  Backscatter  cross  section  for  cylinder  near  plane  as  a  function  of  distance 
AL  at  60°  incidence. 


GammaO 


PC  Cylinder  Above  Conducting  Plane  (70  deg  Inc) 


Figure  13:  Backscatter  cross  section  for  cylinder  near  plane  as  a  function  of  distance 
AZ<  at  70°  incidence. 


PC  Cylinder  Above  Conducting  Plane  (80  deg  Inc) 


Figure  14:  Backscatter  cross  section  for  cylinder  near  plane  as  a  function  of  distance 
A L  at  80°  incidence. 


To  introduce  the  results  of  primary  interest,  we  now  turn  to  example  (c).  A  flat 
surface  alone  scatters  energy  only  in  the  specular  reflection  direction;  however,  if  the 
surface  is  made  slightly  rough,  the  energy  of  the  incident  wave  will  be  scattered  poten¬ 
tially  in  all  upward  directions.  The  scattering  characteristics  of  such  a  rough  surface  are 
shown  in  Figure  15.  Panel  (a)  shows  the  realization  of  a  linear  ocean  surface  at  wind 
speed  U  =  10  m/s.  Panel  ( b )  shows  the  differential  scattering  coefficient  for  a  wave 
(1)  incident  on  this  linear  ocean  surface  alone,  (2)  incident  on  the  linear  ocean  surface 
with  a  cylinder  5  m  below  the  zero-mean  surface  level,  and  (3)  a  PC  cylinder  5  m  below 
a  planar  surface.  The  scattering  characteristics  of  the  linear  ocean  surface  alone  were 
computed  numerically  using  MOM,  which  illustrates  the  technique  of  combining  objects 
whose  individual  scattering  characteristics  are  more  easily  managed. 

Several  points  should  be  noted.  The  effects  of  the  cylinder  interacting  with  the 
rough  surface  exceed  the  rough-surface  level  at  essentially  all  scattering  angles  greater 
than  —20°.  In  the  backscatter  region,  the  scatter  from  the  cylinder  dominates,  but  the 
rough-surface  level  is  several  dB  higher  than  the  exact  computation  for  the  flat  surface. 
It  should  also  be  noted  that  the  bistatic  scatter  mimima  between  40°  and  50°  for  the 
rough  and  flat  surfaces  are  significantly  displaced  from  one  another.  The  actual  surface 
realization  is  shown  in  Figure  15  (a).  One  can  see  that  even  when  the  surface  height  is 
still  a  small  fraction  of  a  wavelength  there  it  exerts  a  significant  effect  on  the  scattering 
characteristics  of  a  nearby  object.  For  small  levels  of  roughness  the  interaction  tends 
to  systematically  increase  the  backscatter  enhancement. 
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(a)  Surface  Height  —  700/5  pts,  U  =  10  m/s 


(b)  Differential  Scatter  Coefficient  —  0INe  =  70  deg 


Figure  15:  (a)  A  realization  of  a  zero-mean  one-dimensional  linear  ocean  surface.  (6) 
Differential  scatter  coefficient  (1)  for  this  surface  alone,  (2)  for  a  PC  cylinder  5  m  below 
this  surface,  and  (3)  for  a  PC  cylinder  5  m  below  a  planar  surface. 


34 


IV  Scattering  from  Near- Surface  Bubble  Clouds 


As  discussed  in  Section  I,  the  bubble  clouds  are  modeled  as  static  impermeable  objects 
with  known  scattering  characteristics.  At  this  level  of  idealization,  there  is  no  serious 
compromise  in  assuming  the  bubbles  to  be  perfect  voids;  moreover,  to  accommodate  the 
surface  interaction  under  the  high  sea  state  conditions  that  cause  the  bubble  clouds,  we 
are  constrained  to  two  dimensions.  Thus,  we  will  perform  our  model  computations  for  an 
impermeable  cylinder  near  an  irregular  surface.  Because  the  low-frequency  cutoff  of  the 
anomalous  scatter  is  likely  to  be  dependent  on  detailed  bubble  characteristics,  we  have 
used  the  single  intermediate  frequency  of  200  Hz  for  all  the  simulations.  We  have  also 
used  only  one  representative  bubble  diameter  of  1  m.  The  variable  parameters  are  wind 
speed,  incidence  angle,  and  the  distance  of  the  bubble  cloud  from  the  surface  height 
reference  level.  As  the  mutual  interaction  equations  are  formulated  herein,  we  must 
have  a  plane  separating  the  surface  from  the  adjacent  scattering  object  representing 
the  bubble  cloud.  Thus,  the  minimum  distance  at  which  we  can  compute  bubble  cloud 
scattering  is  the  depth  of  the  deepest  trough  for  the  particular  surface  realization. 

The  surface  realizations  are  generated  by  applying  the  improved  linear  represen¬ 
tation  method  described  in  [1],  which  allows  us  to  efficiently  convert  linear  Donelan- 
Pierson  ocean-surface  realizations  to  their  nonlinear  counterparts.  The  acoustic  scatter 
simulations  are  formally  equivalent  to  horizontally  polarized  EM  scattering,  but  the 
scattered  signals  are  reflected  from  the  bottom  side  of  the  surface  realizations.  In  [l] 
we  showed  that  the  top-bottom  asymmetry  of  nonlinear  surfaces  can  change  the  cross 
section  by  as  much  as  5  dB  at  high  sea  states.  For  ocean-surface  realizations,  sea  state 
is  the  principal  parameter.  The  simulations  were  run  for  three  wind  conditions:  U  =  10, 
15,  and  20  m/s,  or  20,  30,  and  40  knots.  Although  some  accommodation  must  be  made 
for  the  one- dimensional  surface  realizations,  these  correspond  roughly  to  Beaufort  sea 
states  5,  7,  and  9,  respectively.  Sea  state  9  is  a  strong  gale.  The  reference  surface  bistatic 
scatter  coefficient  is  computed  for  each  wind  speed,  with  typical  values  lying  between 
-50  and  -60  dB-see  Figure  1.  All  the  results  are  presented  as  unitless  differential  scatter 
coefficients.  The  conversion  to  cross  section  units  is  discussed  in  Section  III. 2. 

As  a  reference  for  each  wind  speed  and  incidence  angle  combination,  we  plot  (1) 
the  surface  realization  for  the  first  rough  surface,  (2)  the  bistatic  scatter  coefficient  for 
the  surface  alone,  (3)  the  analytic  bistatic  scatter  coefficient  for  an  isotropic  scatterer 
below  a  plane  surface,  and  (4)  the  bistatic  scatter  coefficient  at  the  distance  near  the 
first  maxima.  The  cylinder  is  located  along  the  z  axis,  with  the  depth  of  the  maximum 
trough  indicated  on  the  reference  plot.  To  complete  each  set,  we  plot  the  backscatter 
coefficient  for  each  realization  as  a  function  of  distance  from  the  common  reference  level 
and  the  corresponding  analytic  result  for  isotropic  scatterers.  The  reference  level  for 
the  rough  surfaces  is  the  constant  height  level  that  results  from  linearly  detrending  the 


surface  segment  as  described  in  [1].  This  reference  level  is  close  to  (within  10%)  but 
not  equal  to  the  mean  surface  height. 

The  computations  are  performed  for  incidence  angles  of  60°,  70°,  and  80°.  For 
positive  incidence  angles,  wave  motion  is  toward  the  observer.  This  mainly  affects  the 
symmetry  of  the  surface  Doppler  spectra,  which  are  not  computed  for  this  preliminary 
static  bubble  cloud  model.  Two  independent  surface  realizations  were  used  for  each 
computation  to  show  the  sensitivity  to  the  detailed  surface  structure.  The  plots  are 
presented  in  order  of  increasing  wind  speed  at  the  end  of  this  section.  The  top-bottom 
asymmetry  of  the  nonlinear  surface  realizations  is  evident  in  the  reference  plots. 

All  the  simulations  presented  here  were  run  at  the  NASA-Ames  NAS  supercomputer 
center  on  a  CRAY  Y-MP  computer.  The  computation  time  is  driven  by  the  two  matrix 
inversions  required.  For  a  given  surface  realization,  a  single  impedance-matrix  inversion 
provides  the  information  needed  to  compute  the  surface  scattering  function  for  all  in¬ 
cidence  angles.  Computation  of  the  nonlinear  surface  realization  requires  only  Fourier 
transformations  and  efficient  vector  manipulations.  Solution  of  the  mutual  interaction 
equations  requires  another  matrix  inversion.  The  primary  check  on  the  integrity  of 
the  computations  is  energy  conservation.  Because  the  surface  and  the  cylinder  are  im¬ 
penetrable  and  lossless,  all  the  acoustic  energy  that  is  directed  toward  the  surface  and 
cylinder  must  be  accounted  for  in  backscatter.  Depending  on  wind  speed,  to  achieve 
better  than  10%  and  ideally  better  than  1%  overall  energy  conservation  required  700 
to  1200  unknowns.  Tapering  the  incident  field  is  also  necessary  to  control  spurious 
edge  effects,  but  as  with  a  real  experiment  the  beam  effects  on  the  cross  section  can 
be  removed.  For  the  data  presentation  here,  however,  we  work  in  relative  rather  than 
absolute  units.  The  individual  data  runs  required  between  8  and  15  minutes  of  CRAY 
Y-MP  time,  although  the  independent  processors  were  not  used  in  parallel.  Parallel 
processing  could  be  used  for  more  efficient  throughput  for  future  Doppler  computations 
where  more  than  30  computations  per  surface  are  needed. 

The  U  =  10  m/s  data  show  the  same  general  features  as  the  flat-surface  isotropic 
model;  moreover,  the  differences  between  the  two  surface  realizations  is  small  within 
10  m  of  the  surface  reference  level,  which  is  the  region  of  interest.  The  main  feature  to 
note  is  that  at  the  location  of  the  deepest  trough  where  the  computation  is  initiated, 
the  backscatter  coefficient  is  generally  near  or  above  the  free  particle  level  of  ~  —  30 
dB.  Thus,  once  the  bubble  cloud  is  at  the  depth  of  the  deepest  wave  trough,  it  can 
cause  a  significant  backscatter  enhancement  over  its  free-surface  value.  The  incidence- 
angle  dependence  acts  to  lower  the  cross  section  at  the  level  of  the  deepest  trough  as 
the  incidence  angle  is  increased.  Large  incidence  angles  are  most  relevant  to  acoustic 
surface  reverberation. 

The  U  =  15  m/s  data  show  significantly  different  details  but  qualitatively  similar 
behavior.  The  60°  data  evidently  are  already  beyond  the  first  maxima  for  a  flat  surface 


at  the  depth  of  the  deepest  trough.  Because  of  this,  the  first  rough  surface  maxima 
occurs  at  10  m  below  the  reference  level.  At  the  higher  incidence  angles,  however,  the 
initial  surface  backscatter  point  occurs  within  the  envelope  of  the  first  maxima  of  the 
analytic  curve.  Once  the  incidence  angle  reaches  70°,  the  principal  characteristics  of 
the  backscatter  enhancement  are  established. 


The  U  =  20*m/s  surface  realizations  have  crest-to-trough  excursions  that  exceed  the 
incident  wavelength,  as  can  be  seen  from  the  surface  realization  in  the  reference  plots. 
Such  large  surface  roughness  provides  a  severe  test  of  both  the  surface-scatter  and 
mutual-interaction  codes.  We  found  that  we  could  maintain  good  energy  conservation 
for  the  computation  of  the  scattering  function  representing  the  surface,  but  in  the 
ensuing  mutual-interaction  computation,  energy  conservation  was  degraded  significantly 
depending  on  individual  surface  realizations.  For  the  U  =  20  m/s  surface  realizations 
shown  in  this  report  these  errors  may  be  as  high  as  23%.  Surprisingly,  however,  these 
errors  were  confined  to  the  60°  and  70°  incidence  angles.  At  80°,  energy  was  conserved 
to  better  than  8%.  The  relative  error  of  the  ratio  of  the  total  scattered  energy  to  the 
incident  energy  is  summarized  in  Table  1. 


Table  1: 

Relative  Energy  Conservation  Error 

U  =  10  m/s 

U  —  15  m/s 

U  =  20  m/s 

• 

60°  incidence 

0.03-0.04 

0.10-0.13 

0.18-0.23 

70°  incidence 

0.01-0.03 

0.05-0.08 

0.12-0.18 

80°  incidence 

0.01-0.03 

0.01-0.05 

0.00-0.08 

®  While  the  U  =  20  m/s  results  are  somewhat  poorer  in  terms  of  overall  energy 

conservation,  we  believe  that  the  main  features  are  correctly  represented.  We  observe 
much  larger  variations  from  realization  to  realization  and  significant  departures  from 
the  simple  lobing  pattern  predicted  by  the  analytic  curves  we  have  used  for  reference.  A 
complete  characterization  of  a  small  object  scattering  near  such  a  rough  surface  clearly 
demands  a  statistical  characterization,  which  is  beyond  the  scope  of  our  current  effort. 
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(a)  Nonlinear  Surface  Height  —  700/5  pts,  U  =  10  m/s 
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Figure  16:  Differential  scattering  coefficients  for  U  =  10  m/s  at  60°  incidence. 


38 


a 


5  10  15  20  25 

Distance  Delta  L  -  m 


Figure  17: 
incidence. 


Backscatter  cross  section  as  a  function  of  distance  for  U 
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GammaO 


Figure  18:  Differential  scattering  coefficients  for  U  =  10  m/s  at  70°  incidence. 
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Donelan-Plerson  Nonlinear  (U=10  m/s,  ©INC=70  deg) 


Figure  19: 
incidence. 


Backscatter  cross  section  as  a  function  of  distance  for  U  —  10  m/s  at  70° 
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GammoO 


(a)  Nonlinear  Surface  Height  —  700/5  pts,  U  =  10  m/s 


Figure  20:  Differential  scattering  coefficients  for  C/  =  10  m/s  at  80°  incidence. 
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Backscatter  cross  section  as  a  function  of  distance  for  U  =  10  m/s  at  80° 
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(a)  Nonlinear  Surface  Height  —  700/5  pts,  U  =  15  m/s 


(b)  Differential  Scatter  Coefficient  —  0IMC  =  60  deg 


Figure  22:  Differential  scattering  coefficients  for  U  -  15  m/s  at  60°  incidence. 


GammaO 


GammaO 


Figure  24.  Differential  scattering  coefficients  for  U  —  15  m/s  at  70°  incidence 


GammaO 


« 


Donelar.-Plerson  Nonlinear  (U=15  m/s,  0IMC=7Q  deg) 


Figure  25:  Backscatter  cross  section  as  a  function  of  distance  for  U  —  15  m/s  at  70° 
incidence. 
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GammaO 


(a)  Nonlinear  Surface  Height  —  700/5  pts,  U  =  15  m/s 


(b)  Differential  Scatter  Coefficient  —  G,*c  =  80  deg 


Figure  26:  Differential  scattering  coefficients  for  i  =  15  m/s  at  80°  incidence. 


Donelan-Plerson  Nonlinear  (U=15  m/s,  Q1MC=80  deg) 


Figure  27:  Backscatter  cross  section  as  a  function  of  distance  for  U  =  15  m/s  at  80° 
incidence. 
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GammaO 


(a)  Nonlinear  Surface  Height  —  700/5  pts,  U  =  20  m/s 


Distance  -  m 


(b)  Differential  Scatter  Coefficient  —  0,NC  =  60  deg 


Figure  28:  Differential  scattering  coefficients  for  U  =  20  m/s  at  60°  incidence. 
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Donelan-Plerson  Nonlinear  (U=20  m/s,  QINC=60  deg) 


Figure  29:  Backscatter  cross  section  as  a  function  of  distance  for  U  =  20  m/s  at  60° 
incidence. 


* 
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(b)  Differential  Scatter  Coefficient  —  0„,c  =  70  deg 
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GammaO 


GammaO 


(a)  Nonlinear  Surface  Height  —  700/5  pts,  U  =  20  m/s 


(b)  Differential  Scatter  Coefficient  —  0,NC  =  80  deg 


Figure  32:  Differential  scattering  coefficients  for  U  —  20  m/s  at  80°  incidence. 
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GammaO 


Figure  32:  Backscatter  cross  section  as  a  function  of  distance  for  U  =  20  m/s  at  80° 
incidence. 
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V  Conclusions  and  Recommendations 


In  Section  III  and  the  appendix  we  developed  the  mutual-interac  ion  method  in  con¬ 
siderable  detail.  The  solution  of  (14)  characterizes  the  multiply- scattered  wave  fields 
between  two  arbitrary  objects  in  terms  of  their  wave-wave  scattering  functions.  The 
equation  has  the  same  form  as  equation  (5.3)  in  [12],  which  was  used  to  calculate  the 
scattering  of  a  dielectric  sphere  above  a  substrate.  If  the  scattering  operator  defined  by 
(13)  is  approximated  by  a  matrix,  (14)  can  be  solved  as  linear  system,  whereupon  the 
scattered  fields  can  be  computed  from  (15),  (16),  and  (17).  These  fields  are  identified 
schematically  in  Figure  2.  To  the  extent  that  the  scattering  functions  are  known,  the 
mutual-interaction  method  provides  a  potentially  exact  formalism  for  computing  the 
scattered  fields.  For  isotropic  scatterers  near  a  plane  reflecting  surface,  the  equations 
admit  algebraic  solutions,  which  are  summarized  by  (35),  (36),  and  (37).  The  degree  to 
which  these  equations  reproduce  the  exact  solutions  for  the  1-m  cylinder  scattering  at 
200  Hz  is  shown  in  Figures  12,  13,  and  14.  The  departures  from  the  exact  curves,  while 
significant,  are  well  within  other  uncertainties  in  the  model.  The  principal  limitation  of 
the  isotropic  flat  surface  model  is  that  it  does  not  accommodate  wind-induced  surface 
roughness. 

To  model  the  scattering  characteristics  of  subsurface  bubble  clouds  in  the  presence 
of  the  highly  irregular  surfaces  that  spawn  them,  we  have  combined  our  previously 
developed  techniques  for  scattering  from  dynamic  nonlinear  ocean  surfaces  with  the 
mutual-interaction  method.  With  the  restriction  to  two  dimensions,  the  simulations 
remain  computationally  intensive  but  within  acceptable  time/ cost  allocations  of  super¬ 
computer  resources.  As  discussed  in  Section  I,  we  have  extracted  the  simplest  possible 
model  of  a  subsurface  bubble  cloud  that  allows  us  to  investigate  its  scattering  character¬ 
istics  in  the  presence  of  a  rough  surface.  We  chose  a  fixed  bubble- cloud  size  independent 
of  wind  speed  and  calculated  its  scattering  characteristics  as  a  function  of  wind  speed, 
incidence  angle,  and  distance  from  the  surface  reference  level.  These  results  are  sum¬ 
marized  in  Section  IV. 

The  simulations  verify  that  a  nearly  impenetrable  meter-sized  object  scattering  in 
the  presence  of  a  rough  surface  can  easily  account  for  the  anomalous  acoustic  surface 
reverberations.  Indeed,  our  backscatter  coefficient  estimates,  relative  to  the  free-surface 
backscatter  coefficient,  are  much  too  high.  This  ratio  is  easily  manipulated,  however, 
by  changing  the  bubble-cloud  cross  section.  Thus,  we  attribute  this  discrepancy  mainly 
to  the  simplicity  of  the  bubble  cloud  model.  The  integrity  of  the  model  depends  more 
critically  on  how  well  it  explains  the  observed  wind-speed  dependence  of  the  anomalous 
reverberations.  In  this  regard,  a  consistent  feature  lias  emerged  from  the  simulation 
summarized  in  Section  IV. 

The  backscatter  coefficient  at  the  level  of  the  deepest  trough,  where  the  scatter 
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Chapman-Harris  200  Hz 


Figure  33:  Comparison  of  Chapman-Harris  curves  (solid)  with  wind-independent  cloud 
model  (dotted). 

calculation  is  initiated,  has  the  essential  characteristics  of  the  Chapman-Harris  curves 
that  summarize  measured  acoustic  surface  reverberations.  We  can  illustrate  this  in  a 
simple  way  by  using  the  flat  surface  isotropic  scatter  model.  Backscatter  strength  in  this 
simple  model  is  controlled  by  a  single  parameter,  namely  h2  in  (35).  The  remaining 
variable  parameters  other  than  wave  number  are  distance  from  the  surface  reference 
level  and  incidence  angle.  The  two-dimensional  model  will  not  reproduce  the  correct 
wavelength  dependence,  but  all  other  aspects  of  the  model  are  generally  consistent  with 
the  three-dimensional  model.  The  incidence  angle  dependence,  for  example,  is  identical 
in  the  two  models.  To  introduce  a  wind-speed  dependence  in  the  flat-surface  equations, 
we  assume  that  the  subsurface  bubble  population  acts  in  aggregate  as  a  large  void 
only  at  the  depth  of  the  deepest  trough.  We  estimate  this  depth  by  using  (1).  With 
this  hypothesis,  the  only  free  parameter  in  the  model  is  scattering  strength.  We  have 
adjusted  this  level  to  match  the  Chapman-Harris  curve  at  a  wind  speed  of  20  knots. 
The  results  are  shown  in  Figure  3'^>SJ 
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The  bubble-cloud  curves  rise  too  steeply,  but  a  constant  bubble-cloud  cross  section 
independent  of  wind  speed  is  certainly  unrealistic.  We  note  also  that  the  incidence- 
angle  dependence  tends  to  saturate  at  the  highest  wind  speed  more  rapidly  than  the 
Chapman-Harris  curves.  This  is  most  likely  due  to  deficiencies  in  the  isotropic  flat- 
surface  model.  We  believe,  however,  that  the  bubble  cloud  plus  rough  surface  model 
has  the  potential  to  fully  explain  the  observed  acoustic  surface  reverberations.  What 
is  needed  is  a  detailed  parametric  study  guided  by  more  information  on  the  observed 
characteristics  of  sub-surface  bubble  clouds,  particularly  their  wind-speed  dependence. 
The  simulations  can  be  performed  efficiently  with  modest  improvements  in  the  scatter 
code.  At  this  point  in  time,  the  restriction  to  two  dimensions  is  not  severe. 

The  numerical  simulations  also  have  the  potential  to  simulate  broad-band  signal 
characteristics.  Once  appropriate  parameters  are  established,  the  computations  re¬ 
quired  are  straightforward.  Here,  however,  Doppler  characteristics  are  important,  and 
good  data  management  will  be  required  in  order  for  computer  resources  to  be  used  effi¬ 
ciently.  In  the  course  of  developing  the  computer  codes  our  initial  effort  has  necessarily 
been  concentrated  overall  on  computational  integrity.  With  some  effort,  the  simulation 
codes  could  be  run  much  more  efficiently. 
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Appendix 


The  Scattering  Function 
of  a  Single  Scatterer 


In  this  appendix,  we  describe  the  method  of  computing  the  scattering  functions  that 
appear  in  the  mutual-interaction  equations  discussed  in  Section  III.  Consider  a  time- 
harmonic  electromagnetic  field  with  temporal  variation  exp (—iut)  propagating  in  an 
unbounded  medium  characterized  by  the  wave  number  k  in  the  presence  of  a  single 
scatterer.  The  scatterer  is  located  at  the  origin  and  is  assumed  to  be  enclosed  in  a 
volume  V,  which  is  finite  in  the  z  direction.  The  incident  wave  field  scatters  from  this 
object,  and  the  total  wave  field  E  can  be  written  as 

E(r)  =  E,(r)  +  E,(r),  (41) 

where  E,  is  the  incident  wave  and  E,  the  scattered  wave.  In  [10]  we  showed  that  the 
scattered  wave  field  is  given  by 

E»(r)  =  k2  //^G(|r-r'|).S(r')dr',  (42) 

where  S  is  the  equivalent  source  function,  and  Q.  is  the  dyadic  outgoing  Green’s  function. 
Hereafter,  dyadics  will  be  distinguished  from  ordinary  vectors  by  an  underline. 

The  source  function  is  related  to  the  total  wave  that  is  impinging  on  the  scatterer, 
which  is  the  incident  wave  E,  for  the  case  under  consideration  because  the  medium 
contains  only  a  single  object.  The  relation  is  given  by 

S(r')=  JJjr  H(r',  r")*E,(r")dr'',  (.43) 

where  H(r',r")  is  a  dyadic  that  depends  on  the  object's  scattering  characteristics.  It 
has  been  shown  that  H  vanishes  identically  for  r'  or  r"  outside  the  object  [13). 
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The  dyadic  Green’s  function  admits  the  two-dimensional  Fourier  decomposition 


‘  III 2^{ 


dK  ( i  exp{ikg(K)\z  —  z'|} 


(2tt)2  l  2 kg(K) 

x  exp{fK-(p  —  p')} 


«(K) 


S(z  -  z') 
k2  ‘ 


where  r  =  p  +  za.z,  az  is  a  unit  vector  along  the  z  axis,  I  is  the  unit  dyad, 
g(K)  = 


[1  -(K/k)2}1^2  K  <  k 
i[(K/k)2  -  1 )1'2  ilK  >  k 

k(K)  =  I  -  Ar2kk, 


and  the  propagation  vector  k  is  decomposed  as 


k  =  K  +  kz(K) a2  =  K  ±  kg(K) a2 


(44) 

(A5) 

(-46) 

(-47) 


with  the  upper  sign  chosen  for  waves  traveling  in  the  positive  z  axis  and  the  lower  sign 
for  waves  traveling  in  the  negative  z  axis.  (For  evanescent  waves,  the  z  direction  of 
propagation  is  defined  as  the  attenuation  direction  to  be  consistent  with  the  analytic 
continuation,  which  will  be  discussed  later.)  For  instance,  /c(K)  in  (-44)  is  evaluated  as 
»e+  =  I  —  A:~2k+k+  when  the  observation  point  r  is  completely  to  the  right  of  the  object, 
i.e,  when  z  >  max{z'}.  The  range  of  integration  in  (A4)  is  over  all  K,  and  our  definition 
of  g(K)  insures  that  each  Fourier  component  in  the  integrand  is  an  evanescent  wave 
that  propagates  in  the  transverse  direction  and  attenuates  in  both  the  ±z  directions. 

The  parameter  g(K )  was  introduced  to  reflect  the  fact  that  we  are  dealing  only  with 
monochromatic  waves  at  a  fixed  frequency  u.  The  free  space  propagation  vector  k  is 
thus  constrained  by  the  relation  |k|  =  k  =  w/c  where  c  is  the  speed  of  light.  Since 
k  can  vary  in  two  dimensions  only,  a  physically  realizable  monochromatic  wave  field 
E(r)  possesses  only  a  standard  two-dimensional  Fourier  transform.  Let  E(K;z)  be  the 
two-dimensional  spectrum  of  this  wave  field  at  the  plane  z  =  const.  It  follows  that 

E(K;  z)  =  E(K;  z0)exp{tkz(K)(z  -  z0)}  =  E(K;  z0)  exp{ikg(K)\z  -  z0 1},  (-48) 

where  Zo  is  a  constant.  It  is  helpful  to  introduce  the  generalized  three-dimensional 
Fourier  transform  for  the  wave  field  with  the  help  of  the  Dirac  delta  function.  Explicitly, 


E(k±)  =  E±(K;  z  =  0)  2tt 6{kz  =f  kg(K)), 


(.49) 


where  the  same  notation  E  denotes  either  the  two-  or  three-dimensional  Fourier  trans¬ 
form  depending  on  the  context  of  its  argument. 
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In  general,  the  three-dimensional  Fourier  decompositions  of  E(r)  and  H(r',r")  can 
be  written  as 


E.(r")  =  ///(2^E*(k'/,)eXp{lk,"-r"}>  (^0) 

-  ///  (^3  ///  (^rjs2(k',  k")  exp{.(k'.r'  -  k".r")}.  Mil) 

In  the  following  we  are  interested  in  the  field  outside  the  scatterer,  hence  the  delta 
function  in  the  integrand  of  (>14)  does  not  contribute.  Substituting  the  above  quantities 
into  (42)  yields 

X  *t(K).H(k,,k")-Ei(kw)  [^(/CJJ^expiitK^p  -  p')} 
x  exp{ikg(K)\z  -  z  |}  exp{i(k'-r'  -  k"-r"  -f-  k'"*r")}.  (.412) 

The  components  of  K,  k',  k",  and  k"'  are  independent  real  variables,  but  the  2  com¬ 
ponent  of  k,  kz,  is  constrained  by  the  relation  k2z  +  K2  =  k2 .  In  addition,  the  ±  signs 
were  conveniently  introduced  in  the  exponent  of  exp{zK-(p  —  p')}  because  G  is  an  even 
function  of  (p  —  p').  Because  the  integrand  vanishes  for  r"  outside  V',  we  can  extend 
the  r"  limits  to  infinity  and  perform  the  integration  over  r"  and  k'"  to  obtain 


-  I  ///.  *■  III.  *•  //  igi  ///  ift  ///  III 51 


x  [g(K)\  1  exp{±fK»(p  —  p')  4-  ikg(K)\z  —  z'\}  exp{fk'-r'}.  (413) 
Upon  resolving  the  absolute  value  sign,  we  obtain 


x[p(X)j  1  exp{t(k'  ±  k(K))-r'}  exp{^tk(K).r}, 


(,114) 


where  the  upper  sign  applies  if  z  <  z'  and  kz  >  0  or  if  2  >  z'  and  kz  <  0  ,  and  the  lower 
sign  applies  if  z  <  z'  and  kz  <  0  or  if  z  >  z'  and  kz  >  0. 

To  perform  the  r'  integration,  we  further  restrict  our  observation  point  to  the  regions 
outside  the  slab  a  <  z'  <  b  where  a  and  b  are  the  minimum  and  maximum  of  the 
projection  of  the  scatterer  onto  the  2  axis.  In  this  case,  the  above  integral  can  be 
rewritten  in  a  compact  form  as 

x  [g(K)\  1  exp{i(k' ±  k(  K  )  )-r'}  e:-:p{-;k(  I\  ur}.  :  a  or  :  b.  (.415) 
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where  now  the  choice  of  sign  depends  on  whether  k(K)  is  directed  away  from  or  toward 
the  origin.  In  the  former  case  the  plus  sign  is  chosen  (outgoing  waves),  and  in  the  latter 
case  the  minus  sign  (incoming  waves).  The  r',  k'  integrations  are  readily  evaluated 
after  the  r'  limits  are  extended  to  infinity  to  yield 

E-(r)  -  f  //  (iy  JJI  (Sp  *(k)-S(±mk),  k")-E.(k“) 

x [g(K)]-1  exp{±fk(K)-r},  z  <  a  ov  z  >  b.  (A 16 ) 


Thus  the  total  scattered  wave  field  at  a  fixed  point  r  is  a  linear  superposition  of 
outgoing  or  incoming  plane  waves.  For  finite  values  of  z,  both  types  of  waves  are 
physically  realizable  for  a  many-scatterer  configuration.  In  the  one-scatterer  problem 
under  consideration  here,  only  outgoing  waves  are  possible  so  that  the  total  scattered 
wave  is  given  by 


Since  Et  is  a  physically  realizable  monochromatic  wave,  the  relation  (.49)  applies  and 
the  above  equation  becomes 


E, 


«  =  7  //  ?5*  //  t£?  *(K)'H(k(K)'  k(K"))-E.CK"; ,  -  0) 


dK 


dK" 


(27r)2  JJ  (2n) 


exp{zk(K)-r} 

P)  ' 


(418) 

This  formula  completely  determines  the  scattered  wave  field  outside  the  scatterer 
once  the  scattering  function  H  is  known.  The  relation  is  simpler  in  the  Fourier  domain. 
In  fact,  taking  the  inverse  two-dimensional  Fourier  transform  yields  the  spectrum  of  the 
scattered  wave  at  a  plane  z  =  const., 


E,(K;  z)  = 


ikex]){ik3,{K)z}  rr  dK" 


2g(K) 


II 


(2tt)- 


!l(K).H(k(K),  k(K"))-E,(K";  z  =  0).  (A19) 


Our  goal,  howeve  ,  is  to  solve  the  integral  equation  (^418)  for  the  scattering  function 
H.  It  can  be  seen  that  the  scattering  function  is  independent  of  the  choice  of  the 
incident  wave.  Thus  we  can  choose  a  unit-amplitude  incident  plane  wave  of  the  form 


E,(r)  —  a,  exp{ik,-r} 

E ,(K";z)  =  at(2r)2S(K"  -  K,)exp{ikt(A',)z}, 


(420) 

(421) 


where  a,  is  the  unit  vector  along  the  direction  of  E,.  and  k,  =  K,  +  a.kz(Kl).  The  K" 
integration  in  (418)  then  yields  a  simpler  integral  equation  for  H. 


(K)'H(k(  K ).  k, )- a, 


exp{?4u(  K  )r}  exp(?K-p) 
g\  A’ ) 


(422) 
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Since  H  is  independent  of  r,  the  above  equation  can  be  solved  by  evaluating  the  integral 
at  large  r  in  a  fixed  direction.  Since  the  condition  z  <  a  or  z  >  b  then  becomes  \z\  — ►  oo, 
we  can  evaluate  this  integral  by  the  method  of  stationary  phase  to  obtain 

ik2  ki)-ai  exp{i(kr  ~  n/A)}/\/8nkr  2-dim 

E,(r)  ~  ^  (>123) 

.  -h2K(Kay H(k„k,)a  t  exp(ikr)/(Airr)  3-dim 

where  k,  denotes  the  scattered  wave  normal  vector  directed  along  r  (i.e.,  k,  =  k^), 
and  K,  is  its  transverse. 

Note  that  the  condition  |z|  — >  oo  excludes  the  case  g(Ks)  7^  0  from  the  above  results. 
Although  the  above  result  was  derived  from  the  argument  of  physical  realizability,  it 
is  valid  for  all  cases  listed  in  (415).  In  fact,  we  can  derive  this  result  independently 
for  each  of  the  four  cases  using  the  method  of  stationary  phase.  Thus,  the  scattering 
function  is  completely  determined  by  the  far-field  behavior  of  the  scattered  wave. 

The  scattering  function  is  also  related  to  the  angular  spectrum  of  the  scattered  waves 
at  a  reference  plane  z  =  z0  outside  the  scatterer.  To  see  this,  we  need  to  evaluate  the 
diffraction  integral 

Es(r)  =  //  ^Ea(K;z0)exp{lkz(K)(z  -  *o)}  exp(zK.p)  (A2A) 

by  the  method  of  stationary  phase  for  large  z.  The  result  is 

kg(K,)E,(K,)exp{i(kr  —  n /A)}/ y/2-rrkr  2-dim 
E,(r)  ~  (A25) 

,  z%(A'J)EJ(K,)exp(zA:r)/(27rr)  3-dim 

where  E,(Kj)  is  a  fictitious  spectral  function  defined  by 

E,(K,)  =  E,(K,;  20)  exp{— ifcz(K,)z0}-  (.426) 

Thus,  we  obtain  the  following  exact  relationship  for  both  two  and  three  dimensions: 

^±(KJ).H(k±(KJ), k,).a,(k,)  =  2<7(tf,)EJ±(KJ)/(zfc),  |K.|  *  k  (427) 

where  for  clarity  the  superscript  ±  has  been  appended  and  the  dependence  of  a,  has 
been  added.  Since  the  above  relation  is  in  the  Fourier  domain,  it  is  understood  that  the 
unit  vector  a,(k,)  characterizes  the  direction  of  polarization  of  one  particular  Fourier 
component  k,. 

Although  E,(K,)  seems  to  be  a  two-dimensional  Fourier  transform  of  the  wave 
field  along  the  plane  z  —  0,  it  is  physically  not.  It  is  obtained  by  propagating  the 
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physical  scattered  wave  E,(K,;2o)  from  the  plane  z  =  z0  back  to  the  plane  z  =  0  as 
if  the  scatterer  were  absent.  Equivalently,  the  scattered  wave  EJ(KJ;20)  appears  as  if 
it  were  generated  from  a  single  scattering  plane  located  at  z  =  0.  Thus,  in  the  two- 
dimensional  Fourier  domain  an  object  can  be  replaced  by  a  transverse  plane  with  the 
same  scattering  characteristics.  For  this  reason  two  physically  different  scatterers  are 
considered  identical  if  they  produce  the  same  spectrum  of  scattered  waves  for  the  same 
incident  wave. 

For  numerical  applications  it  is  convenient  to  introduce  a  vector  version  of  the  scat¬ 
tering  function  defined  by  the  relation 


Wat(K„K,)  =  ^|^iL(K.).fl(k(K.),k(Ki)).ai(k(Ki))> 


|K.|  ^  k.  ( 428 ) 


The  scattering  vector  Sanction  Waj  is  the  two-dimensional  spectrum  of  the  scattered 
waves  produced  by  a  unit-amplitude  incident  plane  wave  polarized  in  the  direction  a,. 
The  scatterer  is  completely  characterized  by  the  scattering  dyadic  H  or,  equivalently, 
three  scattering  vectors  W  associated  with  three  orthogonal  directions. 

In  summary,  we  have  derived  an  expression  for  the  scattering  function  which  is  valid 
for  g(K, )  ^  0.  The  case  where  g{K,)  =  0  can  be  handled  rigorously  by  including  the 
S(z  —  z')  term  of  (.44)  in  the  calculations  from  (.412)  to  (422)  so  that  the  results  will 
be  valid  for  all  z.  We  will,  however,  present  a  short-cut  method  below. 

Determine  the  Scattering  Function  for  Transverse  Propagation 

First,  note  that  H  is  strictly  a  function  of  the  propagation  vectors  k5,  kt,  and  is  therefore 
coordinate-independent.  To  determine  H  for  k,  _L  az,  we  can  rotate  the  coordinates 
(x,y,z)  to  the  new  coordinates  (x,y,  z)  about  the  y  axis  through  an  angle  6.  The 
projection  of  r  onto  the  new  z  axis  is  no  longer  zero,  and  the  stationary  phase  result 
(423)  applies  with  K,  replaced  by  K,,  i.e., 


__  f  Vsirkr  exp{— i(kr  —  37t/ 4) }  2-dim 

«(K,)-H(k4,k,).a, - k~2E,(r)  x  l  ,  (429) 

(  47rr  exp(  — ikr )  3-dim 

which  is  valid  for  k5(Kj)  not  perpendicular  to  the  z  axis.  Since  K,  is  related  to  Ks  by  a 
simple  coordinate  rotation,  the  scattering  function  at  k4(K4)  _L  a;  can  be  obtained  when 
we  appropriately  evaluate  the  above  equation.  Since  the  coordinate  transformation  is 
continuous  in  6  and  all  the  propagation  vectors  are  fixed  in  (429).  we  simply  take  the 
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limit  as  8  — *  0.  This  is  equivalent  to  taking  the  limit  of  (.423)  as  r  — *  rap.  The  result  is 

V8nkr  exp{—i(kr  -  37t/4)}  2-dim 

&.(k*P)‘H{ksip,  k,)*a,  ~  -k~2Es(rap)  x  , 

L  Airr  exp(—  ikr)  3-dim 

(430) 

where  ap  is  a  unit  vector  along  the  transverse  direction. 

It  is  important  to  explore  the  analyticity  of  H  in  the  variable  K,.  The  right  hand 
side  of  (429),  with  K,  set  to  K,,  is  bounded  if  the  scattered  wave  field  falls  off  with 
distance  no  less  fast  than  a  radiation  field.  We  will  assume  this  is  generally  the  case. 
(There  are  exceptions  such  as  scattering  from  an  infinite,  perfectly  conducting  surface.) 
Since  the  above  scattered  wave  field  was  derived  for  an  incident  plane  wave,  which  is 
an  entire  function  in  r,  it  is  analytic  in  any  homogeneous  region  such  as  the  far  zone. 
Rotating  r  through  all  but  transverse  directions,  which  amounts  to  varying  K,,  we  can 
see  that  H  is  analytic  in  Ks  except  on  the  circle  |K,|  =  k.  The  behavior  of  E,(r),  hence 
H  ,  on  this  circle  can  be  inferred  from  (422).  In  fact,  each  point  K,  on  this  circle  is  a 
pole  and  branch  point  of  Es(r).  Equivalently,  the  point  ksz  =  0  is  a  pole  and  branch 
point  when  EJ(r''  is  viewed  as  a  function  of  ksz. 

The  compact  notation  H(k',  k")  introduced  in  (411)  stands  for  H(K',  k'z ,  K",  k") 
where  each  variable  was  independent  and  was  allowed  to  take  any  real  value.  However, 
what  we  have  obtained  in  (423)  are  the  function  H  defined  on  the  shell  |ks|  =  |k,|  =  k 
in  the  k  space  so  that  the  second-  and  the  fourth-slot  variables  are  no  longer  inde¬ 
pendent  variables.  Sometimes,  to  emphasize  this  fact,  we  write  H(k,,kt)  explicitly 
as  H(K,,  fc„(K,),  K,,  A:ti(K,)).  Fortunately  this  special  form  is  the  most  frequently 
encountered  in  our  applications.  For  simplicity,  we  sometimes  drop  the  dependent  vari¬ 
ables  and  write  H(Ka,K,).  However,  this  notation  is  ambiguous  since  to  each  value  K, 
correspond  two  incident  wave  normal  vectors  kf  =  K,  ±  kg[Kt) az  where  the  ±  sign 
refers  to  the  ±z  direction  of  propagation.  Similarly,  there  are  two  scattered  wave  normal 
vectors  k*  =  K,  ±  kg(K,)az.  Thus,  the  short  notation  H±:fc  or  H(Kf,  K*)  can  be  used 
to  unambiguously  denote  the  four  scattering  functions  H(K,,  ±kg(K,),  Kt,  ±kg(Kl)). 
The  scattering  functions  H++  and  H  characterize  the  forward  scattering  behavior  of 
the  object,  and  are  therefore  called  forward  scattering  functions.  The  other  two  scatter¬ 
ing  functions,  and  H~+,  are  appropriately  called  backward  scattering  functions. 

In  summary,  we  can  find  the  scattering  function  of  a  single  object  by  analytically 
or  numerically  solving  the  problem  of  a  unit-amplitude  incident  plane  wave,  exp(ik,-r), 
scattering  from  this  object.  Formula  (423)  or  (427)  then  determines  the  scattering 
function  in  terms  of  the  far-field  behavior  or  the  two-dimensional  spectrum  of  the  scat¬ 
tered  wave  fields.  Repeating  this  procedure  for  the  whole  set  of  k,'s  will  result  in  a 
complete  description  of  M(kj,k,)  for  real  vectors  k5.  k,.  which  can  then  be  resolved 
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into  four  scattering  functions  H±:t. 

Thus  far  we  have  assumed  that  both  k,  and  ks  were  real  in  deriving  the  result  for 
H.  This  implies  that  0  <  KStt  <  k,  i.e.,  both  kaz  and  klz  are  purely  real.  However, 
we  often  need  H  for  all  real  values  of  Ka  and  Kt  to  compute  the  scattered  waves  in 
a  general  problem  where  evanescent  waves  are  involved.  In  the  following  we  will  show 
that  the  relationship  (.,423)  between  the  scattering  function  and  the  far-field  behavior 
of  the  scattered  wave  is  still  valid  for  all  real  values  of  KJit 


Define  the  Scattering  Function  for  Complex  Propagation  Vector  k4 

When  only  K,  exceeds  k,  it  can  be  shown  that  the  derivation  from  ( ^412)— ( ^416),  and 
hence  the  results,  do  not  change  provided  that  we  use  k3Z  =  ±ik[(K,/ k)2  —  l]1'2  =  ih, 
where  h  is  real.  This  can  be  seen  as  an  analytic  continuation  in  the  complex  k,z  plane 
from  the  point  k,z  =  h  to  the  point  k,z  —  ih  along  the  path  hexp(id)  where  6  runs  from 
0  to  7r/2.  This  is  possible  because,  as  previously  discussed,  the  scattering  function  is 
analytic  in  kaz  except  at  the  branch  point  kaz  =  0.  The  path  of  analytic  continuation 
must  not  pass  through  the  branch  point  or  the  branch  cut.  This  branch  cut,  which  is 
implicitly  constrained  by  the  definition  of  g(K )  in  (A5),  can  be  chosen  to  lie  either  in 
the  second  or  the  fourth  quadrant  only. 


Define  the  Scattering  Function  for  Complex  Propagation  Vector  kt 


However,  the  technique  of  analytic  continuation  is  not  obvious  in  the  case  where  kz(Kf) 
becomes  purely  imaginary  since  the  k"  integration  in  (^412)  may  not  collapse  H(k',  k") 
into  H(k',k,).  The  main  difficulty  is  that  analytic  continuation  on  the  variable  k,z 
in  (-A21)  fails  since  the  delta  function  is  singular.  The  delta  function  is  not  strictly 
a  function,  and  any  operation  on  it  must  be  properly  carried  out  on  a  sequence  of 
functions  that  approach  it  in  the  limit.  The  choice  of  this  sequence  is  arbitrary. 


The  following  sequence  and  its  inverse  Fourier  transform  suffice  for  our  purposes: 


Vv {kz;  klz) 

Tpr{z]klz) 


27 rr  1  exp {  — 7rr  2(kz  —  klz)2} 
exp(— r222/47r)  exp (iktzz) 


r  real,  nonzero 


(-431) 


since  in  the  limit  r  — »  0  we  have 


Vv  — ♦  exp(iA:„z),  ifT  — ♦  2ir6(kz  -  klz).  (-432) 

Note  that  both  sequences  are  analytic  in  klz  for  r  at  0.  Thus  they  can  be  analytically 
continued  from  the  real  axis  to  the  imaginary  axis  in  the  complex  kt.  plane: 

V’r (z;ktz  =  ieth)  =  exp( -r2:2  4;r )  exp! -eth:  ) 
rjjT{kz]  kiz  =  ieth)  —  2tt_1  exp{ k.  -  ie,h  )2} 
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where  e,  =  ±1  and  h  =  k\g(Kt)\.  Thus  the  incident  wave  in  ( >120,  421)  can  be 
approximated  by  the  sequence 


E,iT(r)  =  exp(iK-p)  rpT(z\ktz)  (434) 

E,,r(k")  =  (2*)26(K"-Kt)Mk';-,ktz)  (435) 

where  ktz  assumes  a  purely  imaginary  value  ieth. 

Substituting  (>135)  into  (>113)  and  performing  the  k"  integration  we  obtain  a  se¬ 
quence  of  scattered  wave  fields  (cf.  (>114)): 


-  VlL*JI£fUJ 


dK 


(2*)’ 


dk' 

(2tt)3 


x  f  dk"  «.(K)-H(k',  K,,  k")-a.t  r  ^xpf— 7rr  2(k"  —  ieth)2} 

J  —  OO 

x  [g(K)}-1  exp{i(k'  ±  k(K))-r'}  exp{:p2k(K).r}.  (>136) 


The  k"  integration  can  be  evaluated  by  the  method  of  steepest  descent  if  H  is  assumed 
to  be  analytic  over  the  real  and  the  imaginary  axes  of  the  complex  k "  plane.  The 
real  k"  path  of  integration  can  be  deformed  into  a  path  ( C )  which  is  essentially  the 
same  as  the  old  path  except  for  a  portion  near  the  origin  where  it  is  deformed  to 
pass  through  the  saddle  point  ieth.  Most  of  the  contribution  to  the  k"  integral  then 
comes  from  the  neighborhood  of  this  saddle  point.  In  the  limit  of  infinitesimal  r, 
r-1  exy{--KT~2(k"  -  iet/i)2}  becomes  6(k"  -ieth)  and  EJiT(r)  approaches  E,(r)  as  given 
in  (>114).  If  we  repeat  the  calculations  that  follow  Eq.  (.414),  we  will  obtain  the  same 
result  listed  in  (423). 

Note  that  Eq.  (434)  with  r  =  0  represents  an  incident  wave  that  attenuates  along 
the  e,~  direction  but  grows  unbounded  in  the  reverse  direction.  Such  a  wave  does  not 
posses  a  Fourier  transform;  however,  E,  r(r)  with  a  nonvanishing  r  does  have  a  Fourier 
transform. 

Tk  i  analyticity  of  H  in  the  variable  K,  can  be  inferred  directly  from  the  reciprocity 
principle: 

fi(k(Kf),  k(K,))  =  H(-k(K,),  -k(K,)).  (437) 

Consequently,  the  scattering  function  has  the  same  behavior  in  both  variables  K,  and 
K,.  It  is  generally  analytic  in  KJ  t  except  on  the  circle  |K,,,|  =  k.  It  can  be  analytically 
contin  :ed  from  purely  real  values  of  k,Ztlz  into  purely  imaginary  values  ksl  lz  =  ±i(k2  - 
K]  J1'  where  the  choice  of  sign  is  dictated  by  the  attenuation  direction  of  each  wave. 

In  the  following  we  will  study  the  general  scattering  function  for  two-dimensional 
problems  and  consider  specific  scatterers  such  as  cylinders  or  perfect  conducting  planes. 
Some  notations  to  be  used  are  r  =  axx  -  a;c  =  Ir.f))  where  6  is  the  polar  angle  with 
respect  to  the  2  axis  and  k  =  a  xkx  -t-  a  -  k-  —  ( k.O ). 
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Determine  the  Scattering  Function  in  Two  Dimensions 

Without  loss  of  generality,  we  will  show  the  procedure  to  obtain  the  scattering  function 
H(k(A:JX),  k(fctI))  for  all  real  values  of  k3X  and  kxx  for  an  object  located  at  the  origin  in 
two  dimensions. 

Let  us  assume  that  the  problem  of  a  unit-amplitude  incident  plane  wave  with  real 
ki  scattering  from  this  object  can  be  solved  analytically  or  numerically  by  some  method 
such  as  the  method  of  moments.  Let  us  further  assume  that  the  solution  for  the  scattered 
waves  outside  this  two-dimensional  scatterer  can  be  written  as  a  linear  superposition  of 
outgoing  waves  of  the  form 

E.(r)=  £  exp(ma).  (438) 

n=  —  oo 


Here,  an  and  bn  are  the  amplitudes  of  the  nth  normal  mode,  H W  is  the  Hankel  function 
of  the  first  kind  of  order  n,  and  a  =  9,  —  6X  is  the  angle  of  k,  with  respect  to  k,  and  is 
given  by 

a  =  ±  cos-1(k,-kt/A:2)  ( ^439) 

with  the  correct  sign  prescribed  by  a  =  sin_1(k,  x  k ,/k2).  The  scattering  function  can 
be  immediately  inferred  from  relations  (423,  438)  with  the  help  of  the  asymptotic  form 
of  the  Hankel  function  for  large  r  as 


h,x  )•  H( k(  k3X ) ,  k(A:tI))*a, 


_4_ 

ik2 


X!  [aran  +  a  A]  ( -i)n  exp(zna), 

Tl~  — ■  OO 


(440) 


where  —k<  k3XXX  <  k. 

To  derive  H  for  ali  real  values  of  ksx  and  klx  we  need  to  analytically  continue  the 
right-hand  side  of  (440)  from  the  real  axis  to  the  imaginary  axis  of  the  complex  k,z  and 
kXi  planes  along  a  path  in  the  first  and  the  third  quadrants.  The  angle  a  will  become 
complex, 

q  =  Qr  +  iax,  ctj-  and  a,  real,  (^41) 

and  the  constants  an,  bn,  if  dependent  on  k3z(ksx )  and  ktz(klx),  are  assumed  to  be 
analytic  in  these  variables  so  that  they  can  be  analytically  continued.  The  complete 
expression  for  H  can  be  written  formally  as  follows: 


&.(k3X )*il( ksx ,  ksl(ksx'),  kxx ,  kxz(ktx)')'Hl  — 


^  1^2  ^  1  ^x^n( h9X ,  ksz ( k3X ) ,  klx .  klz ( ktx  ) )  ■*-  3zbn{  ksx .  k3Z (  k3X  ) .  klx ,  kiz (  kix  ) J I 


(— i)n exp(— no,) exp(mar ).  — oc  -  k, Xl 


.442 
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Our  goal  is  to  determine  aT  and  a,  for  the  case  where  either  one  or  both  of  k3Z,  klz  are 
purely  imaginary  of  the  form  ktz  =  ie,\kaz\,  klz  =  with  the  sign  e4il  =  ±1  refers 

to  the  direction  of  attenuation  along  the  ±z  axis.  In  general,  we  can  derive  ar  and  a, 
from 


cos  a  =  &~2k,-k,  =  £  +  irj,  <f  and  77  real, 
by  making  use  of  the  identity 

cos-1(£  ±  it])  =  ±|  cos-1s  -  iln(t  +  Vt2  —  1)  j 

where 

t  =  Wti  +  l)2  +  v2  +  $\/ti-i)2  +  v2- 


M43) 

(444) 

(445) 

(446) 


(i)  Compute  a  When  Only  ktz  Is  Purely  Imaginary 
In  this  case  we  find 


• 

(£  ±  l)2  +  772  = 

{ktx  ±  k3X)2/k2, 

(447) 

s  = 

cos{±[0,  -  sgn(fcJX)7r/2]  |, 

(448) 

t  = 

where  sgn(x)  is  the  sign  of  x.  Thus, 

\ktx/k\, 

(449) 

• 

aT  =  cos- 

1  s  =  9,  -  sgn(/t,I)7r/2, 

\klx\  +  \J k2x  -  k2 

(450) 

a,  =  ±  In 

k 

(451) 

Here,  the  plus  sign  in  (.448)  was  chosen  since  it  gives  the  correct  interpretation  of  ar  as 
the  angle  between  the  real  directions  of  propagation  of  the  incident  and  scattered  waves. 
The  sign  of  a,  will  be  determined  on  the  individual  case  basis  for  each  normal  mode  n 
of  the  scattering  function  based  on  whether  the  mode  is  growing  or  is  damped  as  |A;,J 
increases.  Since  H(fcJX,  k,z,  kix,  was  defined  with  the  signs  e,  =  ±1  denoting  the 

±z  direction  of  attenuation,  the  four  scattering  functions  are  always  damped.  Taking 
into  account  the  sign  of  the  normal  mode  number  n,  we  obtain  the  following  results  in 
connection  with  (442): 

aT  =  sgn(k3X)ir  /  2  -  0,, 

exp(inQt)  =  - ~=  -  ■  (452) 

-  V  ^  F 

a0  =  l)  =  b,,. 
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Note  that  the  sign  of  77  in  (.443)  ,  which  is  related  to  the  eaiz  direction  of  attenuation 
of  the  incident  or  scattered  waves,  does  not  play  a  role  in  determining  a  unique  solution 
for  at.  (This  is  also  true  if  we  use  sin  a.)  Had  we  used  the  formula 


a  =  tan  1 


k;  x  k. 


k»-k, 

we  would  have  obtained  the  result 


=  tan 


-1 


|"  foix^sz  "1"  ie%kaxyjk*x  k? 

ktxkax  +  ietkaz^/k^x  -  k2  J 


a,  =  i  In 


|*«  -  e^kl  -  W | 


and  the  natural  choice  of  the  solution  would  have  been 


H 


exp(inai)  = 


■  I ktx  +  e,  \J k2lx  -  k2\. 


,  k  ^  kax  ^  k. 


(453) 


(454) 


(455) 


This  solution  has  great  appeal  since  it  depends  on  the  direction  of  attenuation  (e,)  and  f 

appears  to  be  unique  for  each  normal  mode  n.  However,  it  must  be  discarded  because 

of  its  unbounded  growth  in  one  direction  along  the  ktx  axis.  This  peculiar  behavior 

results  from  the  fact  that  tana,  as  a  function  of  klx  at  a  fixed  value  of  kax  —  0,  has 

two  real  singularities,  ktx  =  ±k  as  seen  from  equation  (453).  Since  the  solution  a  can 

take  on  different  character  on  different  sides  of  the  singularities,  it  is  impropriate  to  use 

formula  (453)  to  derive  a,.  • 


(n)  Compute  a  When  Only  kaz  Is  Purely  Imaginary 


In  this  case  we  obtain  the  same  results  as  listed  in  (447-451)  except  that  the  subscripts 
z,  s  are  exchanged.  The  results  are 


olt 

exp(zna,) 


CL  0 


6,  -  sgn(A:,x)7r/2, 

r _ 

M  k,x\  +  Jk*„-k*/ 

0  =  b0. 


(456) 
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(Hi)  Compute  a  When  Both  klz  And  ksz  Are  Purely  Imaginary 
In  this  case  we  have 


|{±1|  =  sgn(kaxkix)((  ±  1), 
•s  =  sgn  (ksxkix), 


t  =  sgn  {k,xktx) 


k,xkix  ±  Jkl  -  k*Jkl  -  P 


(A51) 

(A58) 

(A59) 


After  some  algebra  we  obtain 


ar  =  sgn(fc«)7r/2  -  sgn(A:tI)7r/2, 
exp(znat)  =  ( - *==)  f - k  == 

\\klx |  +  yjkl  -  k*J  \\ksx\  +  y/kl  -  P 

a  o  =  0  =  6q- 


(A60) 


In  the  following  we  will  apply  the  above  results  to  an  infinite  cylinder,  an  infinite 
conducting  plane,  and  an  infinite  conducting  surface. 

Scattering  Function  of  an  Infinite  Cylinder 

An  infinitely  long  cylinder  of  radius  a  and  refractive  index  m  is  assumed  to  lie  along 
the  y  axis.  The  incident  propagation  vector  k,  is  assumed  to  be  normal  to  the  cylinder 
axis  (k%y  =  0).  We  will  consider  two  cases  of  polarization  separately  and  make  use  of 
the  results  given  in  Chapter  8  of  Bohren  et  al. 

(i)  Axial  Polarization:  For  incident  waves  polarized  along  the  cylinder  axis  we  have  the 
following  results: 


E,||(r)  =  avEt  exp(zkt*r)  =  ayE,  exp(zP  cos  a), 


(-461) 


E.||(r)  ~  r-  exp{z(P  —  37r/4)}  ^  cnexp(zna),  r  — ♦  oo,  (A62) 


Cn  —  C_n  — 


Jn(mka)J^(ka)  —  mJ'n(mka)Jn(ka ) 
Jn(mka)Hn^'(ka )  -  mJ'Jjnka)  Hn\ka) 


,  n  >  0, 


( .463) 


where  a  =  0,  —  Ql  is  the  angle  of  k,  from  kt,  Jn  is  the  Bessel  function  of  order  n,  and 
Hn  is  the  Hankel  function  of  the  first  kind  of  order  n.  Since  the  scattering  function  is 
produced  by  a  unit-amplitude  incident  plane  wave,  we  set  E,  to  unity  and  identify 


*(  k9X ) *£L(  kj ,  k, } ■  Uy  —  au  H\i  (  k$ .  k,  i . 


( .464) 


where,  for  real  values  of  kxz  and  k,z, 

oo  oo 

i/||(kJt  k.)  =  4 ik~2  ^2  cnexp(ina)  =  4 ic0k~2  +  8 ik~2  y~)  cncos(na).  (.465) 

n=  —  oo  n=l 

Since  Cn  does  not  depend  on  kaz  or  kxz,  we  can  analytically  continue  only  a  to  imaginary 
values  of  k,z  or  ktz  and  set  c0  to  zero  to  obtain 

OO 

//y(k4,kt)  =  8 ik~2  cn  exp(mo!1)  cos(mar)  ( >166) 

n=l 

where  ar  and  exp(znctt)  are  given  by  (.452)  if  only  klz  is  purely  imaginary,  ( >456)  if  only 
k3Z  is  purely  imaginary,  or  (.460)  if  both  kxz  and  k,z  are  purely  imaginary. 


( ii )  Transverse  Polarization:  The  incident  wave  is  assumed  to  be  polarized  perpendic¬ 
ular  to  the  cylinder  axis, 

Eix(r)  =  3LtEi  exp(ik,«r)  =  a,E,  exp(ikr  cos  a),  ( >16 7 ) 


where  a,  =  a,  x  k Jk  =  a«t  is  a  unit  vector  pointing  in  the  direction  of  increasing  8X. 
The  results  are  summarized  as  follows: 


E,±(r) 

dn  —  d_n 


/  2  ” 

E'\  nkr  expW*T  ~  ’’V4)}  dn  exp(ina),  r 


mJn(mka)J^(ka )  -  J'n{mka)Jn{ka ) 


mJn(mka)Hn^'(ka )  —  J^(mka)Hnl>(ka) 


(0/ 


n  >  0 


oo  (.468) 
(.469) 


*(*«)*Ii(k(A,I),k(A^)).af(k(A^x))  =  a,(k(A:4X))  H±{k(k3X),  k(A,r))  (>170) 

a,(k(A:tI))  =  av  x  k (ktx)/k  (>171) 

a,(k(fcJX))  =  ay  x  k(k,x)/k  (>172) 

where  H±  is  given  by  the  same  expressions  as  (>165,  >166)  with  cn  replaced  by  dn. 

Note  that  the  unit  vector  a,  or  a§  in  (>167,  >168)  characterizes  the  direction  of  the 
total  incident  or  scattered  wave  at  a  point  r.  However,  the  unit  vector  a,(k(fctx))  or 
a,(k(fcJX))  in  (>170)  is  the  direction  of  each  Fourier  component  of  the  incident  or  scat¬ 
tered  wave.  These  vectors  depend  strictly  on  the  propagation  vector  of  each  Fourier 
component  and  are  defined  by  equations  (.471.  .472)  for  real-  or  complex-valued  vec¬ 
tors  k, ,.  When  |fc4X|lx|  >  k,  the  unit  vectors  become  complex  and  the  corresponding 
evanescent  Fourier  components  become  elliptically  polarized. 
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It  is  important  to  point  out  that  the  cylinder  is  azimuthally  (6)  symmetric  and 
•  thus  the  scattering  function  H  is  expected  to  be  even  in  k3X  for  the  case  k,  =  ±a zk. 

Our  above  results  are  consistent  with  this  expectation  and  confirm  that  the  solution 
unbounded  in  one  direction  of  ksx  derived  from  tana  is  incorrect. 


Scattering  Function  of  an  Infinite  Conducting  Plane 

Consider  a  horizontal,  infinite  conducting  plane  located  at  the  origin  with  its  normal 
directed  along  the  z  direction.  The  incident  wave  is  assumed  to  be  propagating  from 
above  the  surface.  Since  the  tangential  component  of  the  total  wave  field  must  vanish 
on  the  conducting  surface,  we  can  obtain  the  following  solution: 


Ei(r)  =  a.2?,  exp(fk,*r)  =  a,Et  exp{z(K,-p  +  ktzz)} 


(473) 


E  /n  _  |  (-apa,  +  a2a*).at£lexp{i(Kt.p  -  ktzz)}  above  surface 
*  (  — Et(r)  below  surface  ^ 

(Recall  that  the  scattered  wave  is  defined  as  E,  =  Etota{  —  E,.)  Thus  the  total  wave 
below  the  surface  is  identically  zero  as  expected.  Note  that  the  above  solution  is  valid 
in  both  two  and  three  dimensions. 

The  scattering  function  can  be  immediately  found  by  taking  the  two-dimensional 
Fourier  transform  of  the  above  equations  and  setting  Et  to  unity.  The  two  backward 
scattering  functions  are  defined  by  the  solution  above  the  surface  as 


«L+.H+-(K„Kt)  =  «-.H-+(K„K,)  = 


8? T2g(K,) 
ik 


6( K,-K,)(I-2aiai),  (475) 


while  the  two  forward  scattering  functions  are  determined  by  the  solution  below  the 
surface: 


«+-H++(KJ,  Kt)  =  k'.H— (K„K,)  = - 


8t r2g(K3) 
ik 


S(K3  -  K,)I.  (476) 


Scattering  Function  of  an  Infinite  Conducting  Surface 

For  simplicity  let  us  confine  ourselves  to  the  case  where  the  surface  height  varies  only  in 
the  x  direction.  Let  z  =  f(x)  be  the  surface  height  variation  about  the  z  =  0  reference 
level.  Since  electromagnetic  waves  cannot  penetrate  a  perfectly  conducting  surface,  we 
assume  that  the  incident  wave  is  propagating  toward  the  surface  from  above  only,  hence 
only  the  two  scattering  functions  H+“  and  H"“  are  meaningful. 
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Consider  an  arbitrary  downgoing  incident  electromagnetic  wave  which  is  linearly 
polarized  parallel  to  the  surface  (parallel  polarization), 


E«||(r)  =  ayE-{  r), 

or  parallel  to  the  plane  of  incidence  (vertical  polarization): 

E,i(r)  =  a,£-(r), 
a,  =  ay  x  k ,/fc. 


M  77) 


(A  78) 
(-479) 


The  complete  solution  of  the  scattered  waves  for  this  general  problem  has  been  given 
in  our  recent  report  [1].  For  our  purposes  it  suffices  to  summarize  the  results  therein  as 
follows. 

(i)  Horizontal  Polarization:  The  two-dimensional  spectra  of  the  upgoing  (  +  )  and  down¬ 
coming  (  — )  traveling  scattered  waves  are  given  by 

£>+,,.  )  1  f  j  #NJ  ,  ,  IOM, 

Eay(««)  =  ~2  J  Jy(x  > - q(k~) - exp[  —  ik,xx  )dx  ,  (d8U) 


E,Jksx)  =  -Eiy(k,x), 


(-481) 


where  Jy{x)  is  the  effective  surface  current  satisfying  the  integral  equation 

«i(W  -  y  (482) 

By  choosing  a  unit  amplitude  incident  plane  wave,  i?,(r)  =  exp(ik,-r),  we  can  identify 
the  two  scattering  functions 


!S,  (fc»*)*H  {ksx)  kix)’&y  —  nyH^  (k3xxkxx), 

tL  (*,x)-H  (kax,  kxx^-&y  =  ay//,,  (k3X,klx^t 


where 


(-483) 

(-484) 

(-485) 


H+*(k,„k„)  =  MM  (-485) 

A,' ■(*„,*,.)  =  -  M-  (-486) 


(ii)  Vertical  Polarization:  The  results  can  be  summarized  in  our  notation  as  follows: 

EtlJ(  klX  3x  )  =  3i,j(  ) 

®i,j(  ktXi3X )  =  ay  •  k, ( klx  sx  )  ,■ /c 


74 


(-187) 

(.488) 


where 


£.+(M 

£;(*«) 

E~(ktx) 


-5  JJM 


9{k, i)  - 


k3X  df(x') 


k  dx' 

exp{-ikg(k3X)f(x')-ik3Xx'}  , 


-E~(kax), 


fl'(fc.i)  + 


a(k,t) 

ki X  d/(x') 


A;  dx' 


dx',  (>189) 

(>490) 

exp{*fc$(fc„)/(x')  -  iklxx'}  ,  1ft1^ 

- — - ax  .  (>191) 

9\kix) 


Thus  we  can  infer  the  scattering  function  for  this  case  as 

&+(fc.«)-fi+_(*«,  *»W*«)  =  a 3(k3X)H++(k3X,  klx),  (>192) 

^ix )  ~  (^ax)^j_  (^jxj^ix))  (i493) 


where 


Ht+{k3X>  ktx)  =  ^^-E^(k3X), 

Hr(ksx,ktx)  =  2J^E;(k3X)^-^^l6(k3x-klx). 


(.494) 

(-495) 


Conclusion 

In  conclusion,  we  have  derived  the  expression  for  the  scattering  function  in  terms 
of  the  far-held  behevior  or  the  two-dimensional  spectrum  of  the  c-.attered  waves  for  a 
single  object  located  at  the  origin.  This  result  can  generally  be  analytically  continued 
to  imaginary  values  of  the  kz  component  of  the  incident  and  the  scattered  waves  under 
the  assumption  that  the  scattering  function  is  bounded. 
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